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A TWO-DIMENSIONAL 
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By 

James  Sanford  Sirkis 
April,  1988 

Chairman:  Charles  E.  Taylor 

Major  Department:  Aerospace  Engineering,  Mechanics  and 

Engineering  Science 


A two-dimensional  hybrid  experimental-numerical 
technique  for  elastic-plastic  stress  analysis  is  presented. 
This  technique  results  from  merging  two  relatively  new 
technologies  in  engineering  mechanics:  boundary  element 

methods  and  image  processing.  A syntactic  pattern 
recognition  scheme  termed  'Displacement  Pattern  Matching' 
(DPM)  determines  displacement  boundary  conditions  to  be  used 
in  an  elastic-plastic  boundary  element  (EPBEM)  code.  The 
result  is  an  automated  stress  analysis  tool. 
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Displacement  pattern  matching  is  a process  where 
displacements  are  measured  by  tracking  an  arbitrary  array  of 
'black'  spots  on  a 'white'  specimen.  The  digitized  images 
of  the  specimen  are  compared  in  a double  exposure  format  to 
determine  displacements.  Displacement  pattern  matching  is  a 
full  field  technique,  with  spatial  resolution  on  the  order 
of  .0001  inch. 

Displacement  pattern  matching  supplies  the  actual 
specimen  displacement  increments  to  the  von  Mises,  isotropic 
work  hardening,  boundary  element  code.  Given  these 
displacements  and  free  surface  conditions,  EPBEM  is  able  to 
incrementally  calculate  the  internal  state  of  stress  at 
selected  locations.  Results  obtained  for  a variety  of 
geometries  and  loading  conditions  compared  well  with  ANSYS 
finite  element  and  selected  published  experimental  solutions 
and  therefore  are  encouraging. 
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CHAPTER  I 


INTRODUCTION 


Substantiating  postulations  through  experimentation,  or 
using  experimentation  to  draw  new  conclusions,  is 
fundamental  to  science  and  engineering.  This  service  is 
provided  to  theoretical  mechanics  by  the  field  of 
experimental  mechanics.  Before  coming  to  any  conclusions, 
experimentalists  demand  confidence  in  their  work;  as  a 
result,  many  time  consuming  experiments  are  repeated  "just 
to  be  sure."  Experimental  stress  analysis  is  no  different. 
Repeating  experimental  procedures,  often  including  more  than 
one  technique,  is  the  norm  when  searching  for  the  solution 
to  engineering  problems  [1],  The  tedious  nature  of  manual 
data  extraction  and  reduction  led  researchers  to  explore 
automated  extraction  and  reduction  schemes.  Hence,  hybrid 
experimental-numerical  stress  analysis  techniques  evolved. 

The  shear  difference  method  [2]  is  an  early  example  of 
an  experimental-numerical  hybrid  method  in  stress  analysis. 
Since  then,  researchers  have  acknowledged  the  need  for 
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improved  and  more  encompassing  hybrid  techniques  [3].  Prior 
to  the  development  of  digital  image  processing,  most  hybrid 
stress  analysis  techniques  suffered  from  the  need  for  manual 
data  extraction.  Seguci  et  al.  [4]  bridged  the  gap  between 
computer  analysis  and  experimental  data  via  a digitized 
video  picture.  Their  paper  marks  the  beginning  of  the  era 
where  researchers  started  using  digital  image  processing 
techniques  with  structural  analysis  codes  to  produce  viable 
stress  analysis  techniques  [5-6].  During  this  same  era, 
researchers  started  to  exploit  a new  numerical  structural 
analysis  technique,  the  boundary  element  method. 

The  appearance  of  boundary  element  methods  (BEM)  was 
coincident  with  the  advent  of  digital  computers  [7],  but 
they  did  not  become  popular  until  Rizzo  [8]  introduced  the 
direct  boundary  element  method.  Still,  BEM  (sometimes 
called  boundary  integral  equation  methods)  have  yet  to 
achieve  the  acceptance  of  finite  element  methods  (FEM) . The 
reluctance  of  the  technical  community  to  embrace  BEM  is 
largely  attributable  to  its  less  familiar  mathematical 
foundations  [9],  However,  in  the  past  ten  years,  the 
popularity  of  BEM  for  hybrid  and  purely  numerical  studies 
has  greatly  increased. 

Moslehy  and  Ranson  [10],  Umeagukwu  et  al . [11],  and 
Balas  et  al.  [12]  introduced  early  BEM-experimental  hybrid 
techniques,  but  all  require  manual  data  extraction.  Liu 
[13]  introduced  a hybrid  method  using  digital  image 
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processing  and  elastic  BEM.  In  his  method,  digital  image 
correlation  is  used  to  find  displacements  interior  to  the 
boundary.  Liu  then  solves  the  inverse  boundary  element 
problem  to  find  the  boundary  conditions.  Liu's  work  finally 
relieved  the  experimenter  of  data  extraction 
responsibilities . 

Hybrid  elastic  boundary  element-experimental  methods 
are  presently  an  important  research  area  [14-15].  Equally 
important  are  elastic— plastic  boundary  element-experimental 
methods,  but  there  seems  to  be  an  absence  of  research  in 
this  topic.  Elastic-plastic  boundary  element  methods 
(EPBEM)  were  first  conceived  by  Swedlow  and  Cruse  [16]  in 
1971  and  then  implemented  by  Riccardella  [17]  in  1973. 

Since  then,  EPBEM  have  grown  in  popularity.  Telles  [18] 
gives  an  excellent  historical  review  of  the  development  of 
this  numerical  method.  The  rapid  increase  in  the  popularity 
of  elastic-plastic  boundary  element  is  due  to  [18-19]: 

1.  The  BEM  solution  reduces  much  of  the  domain 
by  one  dimension.  That  is,  in  elastic 
regions,  only  the  boundary  is  discretized. 

2.  Only  the  plastic  zone  must  be  domain 
discretized. 

3 . BEM  produces  a reduced  set  of  equations 
(compared  with  FEM) . 

4 . The  user  preparation  to  run  the  codes  is 
greatly  reduced. 

5.  Expertise  requirements  are  greatly  reduced 
from  techniques  dependent  on  mesh  generation. 
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6.  Infinite  and  semi-infinite  problems  can  be 
accurately  and  efficiently  modelled. 

7.  The  displacements,  stresses  and  strains  can 
be  selectively  calculated. 

8.  EPBEM  affords  a high  degree  of  accuracy. 


The  EPBEM  algorithm  used  here  is  conceptually 
eguivalent  to  the  rate  form  of  elastic  BEM  with  a body  force 
[20].  The  governing  differential  equation  for  displacements 
can  be  arranged  so  that  the  homogeneous  part  is  the  Navier 
equation  for  total  displacement  rates.  The  nonhomogeneous 
part  contains  the  nonlinear  terms  associated  with  spatial 
derivatives  of  the  plastic  strain  rate.  These  nonlinear 
terms  constitute  a pseudo  body  force.  The  boundary  element 
solution  approach  leads  to  an  equation  for  the  unknown 
boundary  conditions  in  terms  of  the  known  boundary 
conditions  (assuming  the  body  forces  are  known) . Since  the 
plastic  strain  rates  introduce  an  additional  set  of 
unknowns,  the  boundary  element  method  requires  an  additional 
set  of  equations.  The  constitutive  model  provides  enough 
information  so  that  an  initial  strain  [21]  solution 
procedure  can  be  successful.  This  EPBEM  algorithm  uses  the 
isotropic  hardening,  von  Mises  yield  criterion  discussed  by 
Swedlow  [22]. 


5 


In  the  natural  progression  of  coupling  boundary  element 
methods  and  image  processing  techniques,  this  report 
introduces  a two-dimensional  experimental -EPBEM  technique. 
The  experimental  technique  used  here,  displacement  pattern 
matching,  is  a spin-off  of  the  pattern  mapping  technique, 
first  introduced  by  Fail  [23],  Pattern  mapping  uses  digital 
image  processing  and  syntactic  pattern  recognition 
techniques  to  discern  the  motion  of  spots  applied  to  a 
specimen.  The  history  of  pattern  mapping  can  be  traced 
through  the  photogrid  method  [24-25]  and  the  fine-grid 
method  [26-27].  These  predecessors  measured  strain  by 
measuring  the  local  motion  of  large  grids  of  dots  applied  to 
a model.  Typically,  the  measurements  were  made  by  the 
experimenter  and  thus  subject  to  a large  degree  of  error. 
Kawasaki  et  al.  [28]  automated  the  fine  grid  method  to 
attain  good  accuracy,  but  the  equipment  is  necessarily 
complex  and  dedicated  to  the  fine  grid  task. 


• • • • • 


Figure  1.  Five  by  five  pattern  mapping  grid  with  an 
icon. 
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Pattern  mapping  has  the  advantage  of  using  general 
purpose  image  processing  hardware.  Additionally,  a regular 
grid  need  not  cover  the  entire  specimen.  In  pattern 
mapping,  one  places  an  arbitrarily  shaped  matrix  of  spots 
(for  example  a 5x5  square)  at  the  locations  in  which  the 
displacements  and  strains  are  required  (Figure  1) . Located 
in  this  matrix,  there  is  a special  spot  (termed  an  Icon) 
that  serves  as  a reference  spot  for  determining  the  relative 
locations  of  the  individual  spots  and  the  matrix 
orientation.  The  position  of  each  spot  in  the  matrix  is 
then  recorded  before  and  after  motion.  Differencing  yields 
the  displacement,  and  a strain  definition  (small,  large, 
Lagrangian,  etc.)  provides  a measure  of  strain  [23].  In  the 
proposed  hybrid  stress  analysis,  only  the  displacement  is 
required,  thus  the  name  "displacement  pattern  matching" 

(DPM) . This  requirement  further  relaxes  the  grid  geometry 
restrictions  to  only  placing  a spot  at  the  locations  where 
one  desires  the  displacements.  In  addition,  DPM  does  not 
use  an  Icon,  resulting  in  the  reduction  of  the  syntactic 
recognition  language  from  context  sensitive  to  context  free. 
The  combination  of  displacement  pattern  matching  for  data 
extraction  and  elastic-plastic  boundary  element  methods  for 
structural  analysis  leads  to  a fully  automated  technique  for 
measuring  plastic  stress  and  strain. 
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The  subsequent  chapters  present  several  numerical  and 
experimental  examples.  The  initial  strain  solution 
procedure  is  verified  by  examining  the  expansion  of  an 
elastic-perfectly  plastic  circular  cylinder  loaded  by 
internal  pressure.  This  verification  includes  comparisons 
of  the  EPBEM  solution  with  previously  published  results 
[29].  The  next  investigation  entails  examining  the 
equivalent  stress  concentration  factor  near  the  circular 
cutout  in  an  infinite  plate  under  uniaxial  tension.  This 
sample  problem  verifies  the  elastic,  work  hardening  EPBEM. 

Both  the  elastic-plastic  boundary  element  method  and 
the  displacement  pattern  matching  routines  are  explored 
numerically  in  ways  best  suited  to  reveal  their  character. 
The  influence  of  experimental  error  on  the  boundary  element 
solutions  is  explored  through  parametric  studies.  Rigid 
body  motion  tests  show  the  accuracy  of  DPM,  and  at  the  same 
time,  show  deficiencies  inherent  in  the  numerical  solution 
scheme.  The  numerical  experiments  yield  possible 
explanations  for  the  automated  technique's  system  errors. 

Finally,  three  experiments  explore  the  ability  of  the 
automated  hybrid  technique  as  a nondestructive  stress 
analysis  tool.  All  are  plane  stress  tests  of  1100-H14 
aluminum,  using  a MTS  testing  machine  to  apply  controlled 
loads.  The  first  is  a simple  uniaxial  stress  test  of  a thin 
sheet.  The  second  is  the  investigation  of  a thin  sheet  with 
a central  hole,  subject  to  uniform  end  loads  (perforated 
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strip  tension  test) ; and  the  third  is  a thin  V-notched 
sheet,  subject  to  uniform  end  loads  (notched  tension 
specimen  test) . All  three  experiments  serve  to  show  the 
usefulness  of  this  hybrid  DPM-EPBEM  technique  in  elastic- 
plastic  stress  analysis. 


CHAPTER  II 


THE  ELASTIC-PLASTIC  BOUNDARY  ELEMENT  METHOD 
II. 1 Linear  Elastic  BEM 


The  derivation  of  the  two-dimensional  integral 
equations  governing  linear  elasticity  and  their  solutions 
provide  an  entry  point  to  the  understanding  of  the 
underlying  concepts  of  boundary  element  methods.  The 
elasticity  equations  and  solutions  are  first  derived,  and 
then  they  are  extended  to  include  plastic  flow.  The  basic 
assumptions  in  the  development  are 

1)  The  problem  is  well  posed  (proper 
boundary  conditions  are  specified) . 

2)  The  material  is  isotropic  and  linear 
elastic. 

3)  Small  displacement  theory  is 
applicable. 


Elastic  boundary  element  methods  are  based  on  the 
numerical  solution  of  the  integral  equations  governing  the 
deformation  in  an  elastic  solid.  These  integral  equations 
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can  be  derived  by  starting  with  the  following  form  of  the 
virtual  work  statement: 


where  ti,  u-^,  and  bj[  denote  the  traction,  displacement  and 
body  force  vectors  respectively.  The  stress  tensor  is  a ^ j 
and  the  strain  tensor  is  e ij  . The  differential  dr  is  an 
element  on  the  bounding  curve  and  the  differential  do  is  an 
interior  surface  element.  Also,  the  *-superscript  denotes 
an  arbitrary  equilibrium  set  and  no  superscript  denotes  an 
arbitrary  compatible  set  [30].  The  right  hand  side  of 
Equation  (II-l)  can  be  rewritten  using  an  identical  form  of 
virtual  work  but  with  the  significance  of  the  superscript 
and  no  superscript  interchanged.  The  resulting  equation, 
given  below,  is: 


r 


r 


+ 


bi (P) ui (P) do 


(II-2) 


0 


This  equation  is  known  as  Betti ' s Second  Reciprocal  Theorem 
[18].  In  these  two  equations,  s represents  a point  on  the 
bounding  curve  and  p represents  a point  in  the  interior. 
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Since  the  choice  is  arbitrary,  ui,  t^,  and  bi  are 
chosen  as  the  traction,  displacement,  and  body  force  of  the 
true  body.  The  choice  of  the  *-superscript  set  is  also 
arbitrary;  therefore,  it  is  chosen  as  Kelvin's  solution  to 
the  response  of  an  infinite,  two-dimensional,  elastic  medium 
which  is  subjected  to  arbitrarily  located,  mutually 
orthogonal  point  loads.  Kelvin's  solution  is  found  by  using 
the  Green's  function  approach  of  letting  b*  be  the  two- 
dimensional  Dirac  delta  function  6ijPj(p)  [31];  where  Pj (p) 
for  j =1 , 2 are  two  mutually  orthogonal  unit  loads  (i=l,2). 
This  designation  is  substituted  into  the  adjoint  form  of  the 
Navier  equation  to  yield  the  following  forms  of 
displacement,  traction,  stress  and  strain: 

u*(s)  = Uij (s,p)Pj (p) , 

t*(s)  = T*j (s,p) Pj (p) , 

* , * (H-3) 

eki (P)  = Ejki(s,p) Pj (p) , and 

aki  (P)  = Sjki(s,p)Pj  (p)  ; 

where 

Uj0  (s'p)=  87(1^770  { ( 3-4i/ ) In  (r)  5 i j - r,irfj}. 


-1 


^S,p)  4jt(1-i/)G  ( [ (1-2t/)5ij+2r,ir,  j] 


dr 

dn 


( 1-2 1/ ) [r(  ^nj  -r ^ jni]  } , 
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-1 


8tt  ( 1—  v ) Gr 


{ (l-2i/)  (r^Sij+r,  j5ik) 


(II-4) 


In  Equations  (II-4)  , r2  = (Xp-Xs) 2 +(Yp-Ys)2,  i/  is  Poisson's 

ratio,  ni  is  the  outward  normal  vector  to  the  bounding 

surface,  dr/dn  is  the  normal  derivative,  and  G is  the  shear 

modulus.  The  function  U*j_j(s,p)  is  known  as  Kelvin's 

fundamental  solution;  where  T*ij(s,p),  E*jk-[(s,p),  and 
£ 

S jki (S/P)  are  the  traction,  strain,  and  stress  respectively 
that  are  associated  with  the  fundamental  solution. 
Sokolnikoff  [32]  provides  a detailed  development  of  U*j  and 
most  boundary  element  texts  include  a derivation  of  T* j , 

71*  * 

Ejki/  and  Sjki  (for  example  see  Ref.  [20]). 

Somigliano's  Identity  for  displacements  [33]  is  the 
basis  for  the  boundary  element  method.  This  identity  is 
found  by  substituting  the  fundamental  solution  into  Betti's 
Second  Reciprocal  Theorem. 


ui  (P)  = U*j  (s,p)tj  (s)dr- 


r 


r 


+ U*j  (s)bj  (p)dn. 


(II-5) 


13 


Equation  (II-5)  can  be  restated  as,  given  all  the  boundary 
conditions  (the  body  force  is  always  considered  known) , 
information  at  any  point  inside  the  boundary  is  attainable 
without  finding  the  entire  domain  solution.  In  other  words, 
selective  calculation  of  interior  displacements  is  possible 
by  a line  integral  around  the  bounding  curve.  An 
appropriate  real  estate  analogy  is  that  a buyer  need  never 
see  the  property,  he  only  needs  to  walk  the  perimeter  to 
know  every  detail  of  the  lot  and  the  house.  This  is  a truly 
remarkable  statement.  The  solution  to  Navier's  equation  is 
reduced  in  dimension  (in  the  sense  that  only  the  boundary 
must  be  examined) . The  reduced  dimensionality  (surface  to 
line  integral)  makes  BEM  a very  attractive  numerical 
solution  technique.  But  still,  the  unknown  boundary 
conditions  remain  unknown,  so  the  solution  is  not  yet 
tractable. 

The  above  dilemma  is  easily  remedied  by  taking  the 
limit  of  Equation  (II-5)  as  the  field  point,  p,  approaches 
the  surface  [8],  Mathematically  stated  as 

Lim  u-^(p),  where  St  is  a point  on  r.  (II-6) 

p-»i 

Since  somewhere  on  the  boundary  r=0  (i=s) , and  since 
T*ij(s,p)  is  proportional  to  1/r  and  U*ij  is  proportional  to 
ln(r),  the  limit  in  Equation  (II-6)  produces  a singular 
integral.  This  singular  integral  exists  [8]  and  is 
interpreted  in  the  principal  value  sense  [34],  The 
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existence  proof  of  the  limit  in  Equation  (II-6)  follows  by- 
segmenting  the  integration  path  to  isolate  the  singular 
point,  then  taking  the  limit  of  the  integral  as  the  singular 
integration  path  shrinks  to  zero. 

The  existence  of  the  right  hand  side's  second  integral 
in  Equation  (II-5)  is  verified  by  the  approach  described 
above.  Using  the  notation  in  Figure  2, 


lim 

p-*  0 


T*j  (s,j?)uj  (s)dr=lim 
J r p-  0 


T*j (s,i)uj (s) dr 

J r 


, , * 

+ lim  T^j  (s,  H ) u-;Gdr  ' . 
p~* o . p-r  1 


(II-7) 


The  first  integral  on  the  right  hand  side  of 
Equation  (II-7)  exists  in  the  ordinary  sense  but  the 
integral  around  the  singular  circuit  must  be  viewed 
carefully.  The  singular  integral's  existence  is  confirmed 
by  making  the  following  two  substitutions: 

1)  T*ij  =q(8)/p ; where  g(0)  is  order  one  and  2)  dr=pd0 


lim 

p-*  0 


*2 

g (8 ) de 


02 

g (0) 


dd  . 


(II-8) 


Generally,  the  limit  in  equation  (II-7)  is  not 
evaluated  until  assumptions  about  the  interpolation  of  uj 
over  the  boundary  r are  made  (see  section  I.  3).  Lachat 
[35],  following  closely  the  arguments  of 
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Muskhelishvili  [34],  showed  that  Equation  (II-7)  evaluates 
to  Cj.jUj  regardless  of  the  interpolation  function.  In  the 
strict  sense  C^j  is  a constant  and  is  dependent  on  surface 
smoothness;  but  in  the  numerical  solution,  C^j  is  also 
dependent  on  the  interpolation  function. 

The  same  arguments  hold  for  the  other  two  integrals  on 
the  right  hand  side  of  Equation  (II-5)  . In  the  first 
integral  U*ij  is  proportional  to  ln(p)  and  in  the  third 
integral  the  area  element  is  pdpd<f> . In  both  cases,  the 
limit  evaluates  to  zero.  Assembling  the  aforementioned 
results  yields  the  boundary  element  constraint  equation. 


Cij (s)Uj (s)+ 


Tij (s,i)uj (s) dr= 
r 


Uij (s,i)tj (s) dr 
r 


+ 


(II-9) 


Uij (s,p)bj ( p ) do 

n 
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For  a well-posed  problem,  Equation  (II-9)  provides  the 
unknown  boundary  conditions.  These  boundary  conditions,  in 
turn,  can  be  used  in  Equation  (II-5)  to  give  the 
displacement  at  any  desired  interior  point.  As  an  example, 
if  the  surface  traction  is  prescribed  on  the  surface,  then 
Equation  (II-9)  yields  the  unknown  surface  displacements. 
Then  Equation  (II-5)  gives  the  displacement  solution  in  the 
domain.  This  procedure  is  valid  for  prescribed  tractions, 
prescribed  displacements,  or  mixed  boundary  conditions. 

II. 2 Elastic-Plastic  BEM 


The  governing  equations  in  elastic-plastic  boundary 
element  methods  are  developed  along  the  same  lines  as 
elastic  BEM;  the  main  difference  is  the  idea  of  initial 
strain.  This  section  will  describe  the  initial  strain 
concept  and  its  effect  on  the  BEM  development. 

Additionally,  the  need  for  a second  equation  is  illustrated 
and  a plastic  flow  rule  is  introduced  to  meet  that  need. 

Swedlow  and  Cruse  [16]  were  first  to  formulate  the 
initial  strain  boundary  element  theory.  Since  then, 
Riccardella  [17],  Mendelson  and  Albers  [36],  Mukherjee  [37], 
Telles  [18]  and  others  have  developed  initial  strain 
algorithms  to  solve  engineering  problems.  Although  initial 
strain  theory  is  not  restricted  to  a particular  flow  rule, 
it  will  be  described  in  the  context  of  isotropic  work 
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hardening  and  incremental  plasticity.  In  accordance  with 
accepted  terminology  and  notation,  rate  and  increment  will 
be  used  interchangeably  and  superior  dots  denote  time 
differentiation. 

In  an  elastic-plastic  body,  strains  can  be  categorized 
as  either  elastic  or  plastic.  In  initial  strain  EPBEM,  the 
plastic  strain  is  assumed  to  act  as  a residual  strain;  each 
load  increment  is  just  a deviation  from  that  residual 
strain.  The  adjective  "residual"  refers  to  the  idea  that 
the  solid  knows,  by  some  internal  mechanism,  the  plastic 
strain  increment  before  the  application  of  the  load 
increment.  When  the  increment  is  applied,  the  material  only 
deforms  elastically.  Loosely  interpreted,  the  initial 
strain  concept  states  that  plastic  deformation  occurs  before 
elastic  deformation.  This  idea  of  initial  strain  is 
analogous  to  Martin's  discussion  of  convergence  by  the 
superposition  of  elastic  and  residual  stresses  [38].  As  an 
aside,  initial  strain  has  little  meaning  in  finite  strain 
plasticity  since  numerically,  this  is  a load  and  unload 
process.  Initial  strain,  however,  is  acceptable  in 
incremental  plasticity. 

The  derivation  of  the  integral  equations  governing 
elastic-plastic  BEM  follow  the  same  procedure  outlined  in 
Section  II. 1.  This  time,  however,  the  starting  point  is  a 
rate  form  of  virtual  work; 
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tj[  (s)ui  (s)  dr  + 
r 


bi(p)ui(p)dfi= 

n 


a i j (s  t P) £ i j (P)  dQ  ; (11-10) 

n 


The  superior  dot  denotes  time  differentiation,  the  terms 
designated  by  the  *-superscript  denote  an  arbitrary 
equilibrium  set,  and  the  terms  with  no  superscript  denote  an 
arbitrary  compatible  set.  The  total  strain  tensor  can  be 
decomposed  into  elastic  and  plastic  parts;  thereby, 
rewriting  Equation  (11-10)  as 


* • 

* • 

t^ (s) u^ (s) dr+ 

bi(p)ui(p)dh= 

r 

Q 

* -P 

CTij (s#p) cij (p)dn 

Q 


+ 


(11-11) 


* -e 

CTij  (s,p)  eij  (p)dQ 

Q 


During  each  load  increment  both  elastic  and  plastic 
deformation  are  present.  Even  though  plastic  deformation 
occurs,  the  stress  rate  is  still  related  to  the  elastic 
strain  rate  by  Hooke's  law  [39-40].  Accordingly,  the  stress 
rate  can  be  found  by  isolating  the  elastic  part  of  the  total 
strain  increment  and  using  Hooke's  Law: 

Sij  = 2G (« 1 j - l?j)  + f^(?kk  (11-12) 


Equation  (11-12)  allows  the  following  substitution  to  be 
made  in  equation  (11-11) : ^ij^ij  = «ij^ij. 

Equation  (II-ll)  can  then  be  written  in  a form  similar  to 
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Equation  (II-2)  by  applying  the  rate  form  of  virtual  work 

(with  the  significance  of  the  superscripts  interchanged) 

* • 

to  the  eijo^ij  term; 


t*(s)ui(s)dr 

r 


ti(s)u*(s)dr 

r 


+ 


b*(p)ui(p)dfi 

n 


(11-13) 


bi(p)ui(p)dQ  = 

n 


* • P 

CTij (s,p) efj (p)dn. 

n 


The  difference  between  this  integral  equation  and 
Equation  (II-2)  is  the  integral  which  includes  the  plastic 
strain  rate  term.  As  in  Section  II. 1,  the  terms  with  the 
*-superscript  are  chosen  as  Kelvin's  solution 
(Equation  (II-3) ) and  the  terms  with  no  superscript  are 
chosen  as  the  displacement  rate,  traction  rate,  body  force 
rate  and  strain  rate  of  the  true  body.  These  designations 
lead  to  Somligliano ' s Identity  for  inelastic  materials; 


uj(s)  = 


Uij (s,p)ti(s)dr  - 


Tij (s,p)ui(s)dr  + 


Uij (s,p)bi(p)dh  + 


(11-14) 


* .P 

CTij  (s,p)  efj  (p)dn 


Q 


The  limiting  process  to  find  the  elastic-plastic  boundary 
element  constraint  equation  is  identical  to  that  described 
in  Section  II. 1.  The  plastic  strain  rate  integral  exists 
and  has  no  contribution  to  Cij . 


20 


cij(s)  Uj(s)  + 


Tij (s,i)ui(i)dr  = 


ft 


Uij(s,i)  ti(i)dr 


ft 


(11-15) 


Uij  (s,  Z ) k>i  (i  ) dft  + 


S j ki  ( s / & ) £ ^ ( -0 ) dft  ; 


where  a*ki(s,p)  = S*jki (s,p) Pj (p)  is  used. 

If  the  plastic  strain  is  known  apriori  and  if  the 
problem  is  well  posed,  then  Equation  (11-15)  yields  the 
unknown  boundary  conditions.  The  boundary  information  could 
then  be  used  to  find  the  displacements  in  an  elastic-plastic 
solid  (Equation  (11-14)).  The  problem  here  is  that  prior 
knowledge  of  the  plastic  strain  rate  is  just  an  idealization 
of  the  initial  strain  solution  approach.  In  reality,  the 
plastic  strain  rate  is  an  unknown  quantity;  therefore,  an 
additional  relation  is  needed. 


II. 3 Elastic-Plastic  Flow  Rule 


An  isotropic  work  hardening,  von  Mises  flow  rule  is 
introduced  to  satisfy  the  above  requirement  for  an 
additional  relation.  Although  a multitude  of  flow  rules 
have  been  proposed  [43],  an  isotropic  work  hardening  flow 
rule  provides  an  easily  understood  constitutive  model.  It 
also  provides  a reasonable  approximation  to  material 
response  while  loading  takes  place.  However,  one  should  be 
aware  that  the  von  Mises  yield  criterion  is  limited  in  that 
it  does  not  provide  for  hysteresis  [44]. 
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Figure  3.  Isotropic  hardening  yield  surface. 

The  yield  function  of  an  isotropic  hardening  material 
mathematically  describes  a yield  surface  that  maintains  its 
shape  but  increases  in  size  as  a function  of  the  plastic 
deformation.  As  an  example  (Figure  3),  f(g)  defines  the 
yield  surface,  where  g is  some  measure  of  plastic 
deformation.  The  loading  path  AB  is  purely  elastic,  where  B 
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is  an  infinitesimal  amount  inside  the  yield  surface.  As  the 
solid  is  loaded  along  BC,  plastic  deformation  commences  and 
the  yield  surface  expands.  Now,  unloading  along  CD  is 
purely  elastic  until  D is  reached.  As  a result,  the  elastic 
domain  is  enlarged. 

There  are  infinitely  many  forms  of  isotropic  work 
hardening  to  choose  from;  Koiter  [45]  presents  a general 
development  of  work  hardening,  Hill  [46]  presents  a detailed 
discussion  of  work  hardening  and  its  applications,  and 
Malvern  [47]  gives  a concise  review  of  the  prevalent  models. 
The  isotropic  work  hardening  rule  used  here  employs  a von 
Mises  yield  criterion;  where  the  yield  function  is  single 
valued  and  independent  of  hydrostatic  pressure.  The  yield 
criterion  may  be  stated  as 


f (*ij)-F(J2)>0; 


(11-16) 


where  f is  the  yield  function  and  J2  is  the  second  invariant 
of  the  stress  deviator. 

Stability  in  Drucker's  sense  leads  to  the  statement  of 
normality  [48]  which,  in  turn,  leads  to  a functional  form  of 
the  proposed  constitutive  model.  Normality  is  expressed  as 


A 


' dt 

acTij 


(11-17) 


where  A1  is  a proportionality  factor  that  depends  on 
J2  and  dJ2  [41-42].  Once  on  the  yield  surface, 
f(aij)  = F(J2 )/  additional  plastic  deformation  requires 
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f (aij ) >f(j2) • Therefore  one  may  write 


Aaf_ 
daij  3ffmn 


amn' 


(11-18) 


where  A is  a different  proportionality  factor  and  is  a 
function  of  J2  [49].  Normality  produces  a linear  relation 
between  plastic  strain  rate  and  the  stress  rate.  In  order 
to  make  Equation  (11-18)  useful  for  computations,  specific 
assumptions  need  to  be  made  about  the  form  of  F(J2) . Under 
the  assumption  of  von  Mises  yielding,  the  yield  function  is 
the  von  Mises  equivalent  stress  (aeq)  and  is  defined  as  [47] 


req  - sij  sij 


(11-19) 


where  Sij  = aij  - Sijoj^  is  the  stress  deviator.  The  3/2 
factor  in  Equation  (11-19)  ensures  ae f equals  the  yield 
strength  when  uniaxial  load  is  applied.  For  a work 
hardening,  von  Mises  material,  an  equivalent  plastic  strain 
rate  can  be  defined; 


<eq  = J § eP.  € P j (11-20) 

such  that  dWp  = aeqeeg  where  dWp  is  the  plastic  work 
increment  [44].  The  multiplication  of  Equation  (11-18)  by 
aij  and  substitution  of  dWp  = aeqeeq  [44]  into  the  result, 
leads  to  the  following  expression  for  A. 


a e 
eg  eg 


If 

d a 


mn 


'mn 


af 

d a 


ID 


aij 


A 


(11-21) 
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. at  2 

Then,  by  noting  — damn  = dJ2  and  3J2  = a eq  the 

a crmn 

proportionality  factor  is  simplified  to 


A 


9 


(11-22) 


eq 


where  E^  is  the  tangent  modulus  [50]  of  the  equivalent 
stress-equivalent  strain  curve.  The  definitions  of 
equivalent  stress  and  strain  are  such  that  under  uniaxial 
loading  E>p  is  also  the  tangent  modulus  of  the  uniaxial 
stress-strain  curve. 

Finally,  substituting  Equation  (11-22)  into  (11-18) 
leads  to  a specific  form  of  the  isotropic  flow  rule,  which 
gives  a relation  between  the  plastic  strain  rate,  stress 
deviator  and  the  stress  rate. 

!p  9 sij  smn  z 

• • = / T T _ O 


Elastic  deformation  is  included  by  combining 

Eguation  (11-23),  Hooke's  law,  and  the  total  strain  rate 
e p 

tensor  (de^j  + de^j)  to  yield  a relation  for  the  total 
strain  rate 
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and  the  stress  rate 

3G  S. . Smn  ?mn 

] - — — u ; (H-25) 

o\f  (1+ET/3G) 

where  incompressibility  (ek]c=0  ) is  used.  By  substituting 
Equation  (11-25)  into  Equation  (11-23)  the  final  form  of  the 
von  Mises,  isotropic,  work  hardening  constitutive  model  is 
obtained. 


il 


= 2G[e 


ID 


ID 


?£k 


•P 

6 i j 


3 sij  smn  *mn 
2 a2  (1+  Et/3G) 


(11-26) 


II. 4 Numerical  Implementation  of  EPBEM 

Simple  domains,  simple  boundary  conditions  and 
knowledge  of  the  loading  history  are  all  needed  if  exact 
solutions  to  the  elastic-plastic  boundary  element  equations 
(Equations  11-14,  11-15  and  11-26)  are  to  be  found.  Some 
exact  solutions  are  documented  [51],  but  seldom  are 
researchers  so  fortunate. 

This  section  will  discuss  an  iterative  solution  to 
the  elastic-plastic  boundary  element  problems.  In  doing  so, 
the  surface  and  domain  discretization,  the  interpolation 
functions  used  to  approximate  the  traction  rate  (tjj  , 
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displacement  (ujj  rate,  and  the  plastic  strain  rate  (eP-^j) 
are  presented.  The  system  matrix  formulation  is  discussed 
in  detail  which  includes  the  calculation  of  diagonal 
elements,  the  calculation  of  surface  stresses,  and  the 
handling  of  geometric  discontinuities.  Finally,  the 
iterative  solution  scheme  for  the  elastic-plastic  problem  is 
detailed.  In  order  to  simplify  the  following  discussion, 
body  forces  are  neglected.  Also,  for  convenient 
referencing,  Equations  (11-14),  (11-15),  and  (11-26)  are 

rewritten  below. 


uj  (P) 


U*j (s,p) ti(s)dr  - T*j (s,p)ui(s)dr  + 


r 


(11-27. A) 


cij(s)  Ui(s)  + 


r 


+ Sjki(s,J’)|P.(i)dh. 


r 


(11-27. B) 


eg 


(1+  Et/3G) 


(11-27. C) 
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The  daily  advances  made  in  memory  and  speed  of  modern 
computers  amaze  even  the  most  avid  computer  enthusiasts. 
However,  today's  computers  have  nowhere  near  the  capacity  to 
exactly  model  the  functions  and  domains  commonly  encountered 
in  many  mechanics  problems.  For  this  reason,  domain 
discretization  and  function  interpolation  have  become  a 
science  unto  themselves.  Boundary  element  methods,  just 
like  finite  element  methods  (FEM)  and  finite  difference 
methods  (FDM) , need  a certain  amount  of  surface  and  domain 
quantization.  Fortunately,  the  integral  approach 
significantly  reduces  the  amount  of  discretization 
necessary.  As  seen  in  Equations  (11-27),  both  surface  and 
domain  integrals  must  be  evaluated.  But  in  BEM,  the  domain 
integrals  are  only  evaluated  in  the  plastic  zone  and 
therefore  represent  a significant  advantage  over  FEM  and 
FDM.  Another  important  difference  is  that  the  domain 
discretization  is  used  only  to  evaluate  the  domain  integral. 

A further  difference  between  finite  element  methods  (finite 
difference)  and  boundary  element  methods  lies  in  the 
solution  approach.  FEM  and  FDM  approximate  the  solution 
over  the  domain,  so  the  solution  is  highly  dependent  on  the 
quantization  and  function  interpolation.  In  contrast,  the 
boundary  element  governing  equation  is  solved  exactly  and 
all  approximations  are  made  in  the  boundary  conditions. 
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Figure  4.  Boundary  modelled  with  straight  line 
segments . 

Some  cases  have  been  reported  where  BEM  has  afforded  as 
much  as  an  order  of  magnitude  improvement  in  accuracy 
over  FEM  [52]. 

The  boundary  element  method  requires  some  way  of 
describing  the  bounding  curve  of  a region;  here,  boundary 
discretization  is  achieved  with  straight  line  segments 
(Figure  4) . Straight  line  segments  can  simply  and 
accurately  approximate  a boundary;  in  many  cases  (plates, 
beams,  etc.)  the  approximation  is  exact.  If  necessary,  the 
surface  can  be  more  closely  modeled  by  increasing  the  number 
of  segments.  It  should  be  noted  that  there  is  a practical 
limit  to  the  size  and  number  of  these  straight  line 
segments,  but  that  limit  is  rarely  encountered. 
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Now  assumptions  must  be  made  about  the  behavior  of  the 
traction  and  displacements  over  these  straight  line 
segments.  By  definition,  the  tractions  are  derived  from  the 
spatial  derivatives  of  displacement.  Although  consistency 
in  derivatives  is  often  ignored,  Brebbia  and  Walker  [53] 
postulated  and  Liu  et  al.  [6]  later  showed  that  a consistent 
boundary  element  solution  greatly  exceeds 


Figure  5.  Representative  boundary  element. 

comparable  inconsistent  routines  in  both  computational  speed 
and  accuracy.  In  addition,  when  only  displacement  boundary 
conditions  are  prescribed,  upwards  to  a 50%  increase  in 
speed  can  be  achieved  [6].  This  last  attribute  is  of 
particular  interest  in  hybrid  experimental-numerical 
methods . 

The  proposed  EPBEM  algorithm  strives  for  a balance 
between  simplicity  and  accuracy.  To  this  purpose,  the 
boundary  displacements  are  assumed  to  vary  linearly  along 
the  boundary  segments.  Correspondingly,  the  tractions  are 
assumed  constant  over  each  segment.  In  BEM  jargon  these  are 
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called  constant  (traction)  and  linear  (displacement) 
elements.  Specifically,  the  boundary  is  broken  up  into  N 
elements,  each  element  is  joined  to  its  neighboring  elements 
at  nodes.  The  linear  variations  between  nodes  can  be 
described  by  (Figure  5) 

Ui(x)=Ui(l+2s/L)+Ui(l-2s/L)  for  -L/2  < s < L/2 ; (11-28) 

where  the  superscripts  1 and  2 denote  node  1 and  2 
respectively,  and  the  subscript  i is  1 for  x-displacement 
and  2 for  y-displacement . The  form  of  Equation  (11-28) 
dictates  that  the  displacement  information  is  located  at  the 
nodes,  whereas  the  constant  element  traction  information  is 
located  at  the  mid  points  of  each  segment.  As  an  aside, 
more  sophisticated  discretization  schemes  are  available 
[13,54] . 

The  plastic  strain  rate  still  must  be  quantized 
before  the  numerical  procedure  can  be  described.  The 
plastic  zone  is  discretized  much  like  FEM.  Triangular  cells 
are  chosen  with  the  plastic  strains  assumed  to  be  constant 
over  each  cell.  This  may  seem  identical  to  FEM  but  in 
contrast  to  FEM,  the  cells  are  used  only  to  evaluate  the 
domain  integral.  If  the  interpolating  functions  are  given, 
then  the  constraint  equation  (11-27. B)  and  interior 
displacement  equation  (11-27. A)  can  be  rewritten  as 
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(11-29. B) 


Equation  (11-29. A)  represents  a system  of  2N  equations 
for  2N  unknowns.  Its  implementation  is  best  described 
through  Figure  6.  Each  set  of  two  rows  produced  by  Equation 
(11-29. B)  is  a result  of  placing  the  source  point  (xs,  ys) 
at  the  center  of  an  element  and  then  integrating  its 
influence  around  the  surface.  The  source  point  is  then 
moved  from  element  to  element  until  the  circuit  is  complete. 
Two  rows  are  produced  because  Equation  (11-29. B)  represents 
both  x and  y directions.  The  above  procedure  leads  to  a 
system  of  equations  in  the  form  of 


[T] (u}=[U] (t}+[S] {eP} ; 


(11-30) 
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Figure  6.  Depiction  of  BEM  solution  methodology. 

where  [T]  is  the  displacement  coefficient  matrix,  [U]  is  the 
traction  coefficient  matrix,  (u)  are  the  boundary 
displacement  rates  (some  of  which  are  known) , and  {t}  are 
the  boundary  traction  rates  (some  of  which  are  also  known) . 
The  plastic  strain  coefficient  matrix  is  [S],  and  (eP) 
contains  the  plastic  strain  rates  (all  of  which  are 
unknown) . Grouping  the  known  boundary  conditions  and 
plastic  strain  rate  information  on  the  right  hand  side 
reduces  Equation  (11-29)  to 

[A] (v}={q}+[S] (Ip) ; 


(11-31) 
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where  [A]  is  the  system  matrix,  {v}  contains  the  unknown 
boundary  conditions,  and  {q}  contains  the  known  boundary 
conditions  multiplied  by  the  appropriate  coefficients.  As 
will  be  discussed  later.  System  (11-31)  can  now  be  solved 
for  the  unknown  boundary  conditions. 

There  is  little  difficulty  in  using  a Gaussian 
numerical  quadrature  [55]  to  evaluate  most  of  the  line 
integrals.  However,  care  must  be  exercised  when  evaluating 
the  elements  containing  the  singularity  (diagonal  elements) , 
as  well  as  evaluating  the  domain  integrals.  The  diagonal 
elements  exist  but  are  singular;  therefore,  either 
analytical  or  special  numerical  techniques  are  necessary. 

The  domain  integral  takes  special  care  since  there  is  a 
certain  amount  of  controversy  over  the  derivative  of 
Equation  (11-29. B).  Each  integral  will  be  discussed  in 
detail  since  understanding  the  reasons  they  merit  special 
treatment  is  vital  to  evaluating  the  integrals  accurately. 

A popular  method  for  evaluating  the  diagonal  elements 
will  be  explained  to  emphasize  why  this  method  is 
inappropriate  for  use  in  elastic-plastic  boundary  element 
methods.  When  the  material  is  purely  elastic, 

Equation  (11-30)  becomes 

[T] (u)=[U] { t) . (11-32) 

Cruse  [56]  realized  that  artificial  rigid  body  motion  would 
not  change  the  physical  nature  of  a mechanics  problem  but 
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would  give  the  diagonal  elements  directly.  By  artificially 
applying  a traction  {t}  in  such  a way  as  to  cause  unit  rigid 
body  displacement  {I},  then  removing  the  traction,  the  net 
result  is 

{ t}=0  and  {u}={I}.  (11-33) 

Since  the  off  diagonal  elements  are  numerically  integrable, 
Equation  (11-33)  can  be  solved  for  the  singular  elements. 


2N 

Tii  = "2  Tij  ; (11-34) 

j=l 

where  Tii  incorporates  both  the  singular  integral  and  C^j . 
Had  there  been  a body  force  vector,  say  (b), 

Equation  (11-33)  would  appear  as 

[T] { I+d)={b) . (11-35) 

The  application  of  unit  displacements  would  eliminate  the 
tractions  but  the  additional  displacement,  (d),  from  the 
body  force  is  still  present.  Since  (d)  is  unknown, 

Equation  (11-34)  does  not  help  in  finding  the  diagonal 
(singular)  elements. 

The  EPBEM  constraint  equation  numerically  acts  like  an 
elastic  problem  with  a body  force;  therefore,  the  rigid  body 
motion  technique  is  not  appropriate.  This  leaves  two 
practical  alternatives.  The  first  alternative  is  to 
surround  the  singular  point  with  many  small  sub-elements  and 
then  integrate  numerically  [57].  The  second,  and  preferred, 
method  is  to  integrate  analytically  along  each  singular 
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segment  [17,36].  In  any  case,  the  principal  value 
coefficient  matrix  (Cij)  is  calculated  analytically  [17] 
(shown  below)  and  added  to  the  diagonal  terms  to  form  the 
system  matrix. 


C11  = 


a _ COS  ( 2-1  ] sin  (a  ) 

2n  An  ( 1—w ) 


sin  ( 2-y ) sin  (a ) 
C12  = C21  = An  (l-i/) 


c22  = 


a _ cos  ( 2-y  ) sin  (a  ) 
2 n An ( 1—u ) 


(11-36) 


Figure  7 shows  that  a is  the  corner  angle  and  y is  the  angle 
between  the  bisector  of  a and  the  x-axis. 


Enlargement  of  a corner  used  to 
calculate  Cij . 


Figure  7. 
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The  domain  integration  scheme  used  in  this  algorithm  is 
chosen  to  avoid  the  current  controversy  over  the  correct 
expression  for  the  spatial  derivative  of  Equation  (11-27. A). 
This  derivative  is  a dominating  aspect  since  it  defines  the 
strain  within  the  body.  The  ongoing  debate  revolves  around 
the  existence  and  the  proper  definition  of  the  derivative  of 
the  plastic  strain  rate  integral.  Although  well  documented 
by  Telles  [18]  in  his  introductory  chapter,  the  salient 
features  of  the  controversy  are  presented. 

The  strain  at  any  point  in  the  domain  is  defined  by 

• • • 

2 e i j = (uifj  + u j f jj  . Many  early  researchers  failed  to 
recognize  that  the  integration  path  around  the  domain 
singularity  changes  with  the  load;  this  path  change 
prevents  pulling  the  derivative  directly  under  the  plastic 
strain  rate  integral  [58-61].  Bui  [62]  was  the  first  to 
recognize  this  error  and  formulated  the  correct  derivative; 
Telles  and  Brebbia  [63]  expanded  Bui's  ideas  and  then  later 
carried  them  out  [21] . The  above  techniques  are  both 
elegant  and  independent  of  the  plastic  strain  rate  shape 
function  but  they  are  sometimes  difficult  to  use. 
Specifically,  high  aspect  ratio  cells  cause  the  integration 
scheme  to  give  erroneous  results.  Since  the  EPBEM  algorithm 
presented  here  is  of  the  constant  cell  type,  the  general 
approach  championed  by  Telles  and  co-workers  will  be 
foregone  in  lieu  of  an  implementationally  easier  technique. 
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Figure  8.  Scheme  used  in  domain  integration. 

The  assumption  that  plastic  strain  rate  is  constant 
over  each  interior  cell  is  not  only  computationally  simple, 
it  also  makes  the  above  discussion  inconsequential.  The 
constant  cell  approximation  enables  direct  evaluation  of  the 
plastic  strain  rate  integral  [17,36,37]  by  subdividing  the 
cell  (Figure  8) . The  integrals  over  triangle  1 and  triangle 
2 are  algebraically  added  to  the  integral  over  triangle  3 to 
give  the  desired  integral  over  the  cell  (triangle  4) . 


38 


The  unit  directions  are  defined  so  that  the  algebraic  signs 
are  correct.  Since  the  plastic  strain  rate  kernal  is 
singular  at  the  base  point,  a small  region  surrounding  this 
point  is  excluded  in  the  kernal  evaluation  (Appendix) . This 
integration  scheme  is  equally  valid  when  the  source  point 
coincides  with  the  load  point  (Figure  9) . No  great 
difficulty  exists  in  evaluating  the  integrals;  and  the 
question  surrounding  the  existence  of  the  derivatives  is 
moot  since  plastic  strain  rates  are  found  via  direct 
differentiation. 


Figure  9.  Singularity  inside  the  integration 
cell . 

Another  point  of  computational  interest  is  the 
evaluation  of  stress  rates  along  the  boundary.  Since  Sj]^ 
has  a singularity  0(r“2),  the  integral  approach  to  calculate 
stress  rates  near  the  boundary  introduces  large  errors 
(since  r«l  near  the  boundary) . Instead,  the  EPBEM  employs 
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the  technique  sketched  by  Hartmann  [64].  Although 
introduced  in  the  context  of  elastostatics , the  method  is 
equally  valid  for  elasto-plasticity . There  are  seven 
unknown  quantities  on  any  smooth  boundary  point  in  a two- 
dimensional  body,  three  stress  rates  j and  four 
derivatives  Uj^j.  There  are  also  seven  equations  available: 

1 ) Hooke ' s Law 

*ij  = 2G?ij  + ?kk  (11-37. A) 

2)  The  definition  of  traction  rate 

j = ^ i j nj  (11-37  . B) 


3)  The  derivative  of  along  the  arclength. 

du^  . dx . 

= Ui,X-i  — 1 

ds^  1 J ds 2 


(11-37. C) 


Note  that  Equation  (11-37. A)  uses  only  the  elastic  portion 
of  the  strain  rates  so  as  to  keep  the  number  of  unknowns  at 
seven;  Equation  ( 11-37. C)  uses  only  the  elastic 
displacements.  Since  displacements  are  approximated  by 
linear  functions,  the  derivatives  along  the  arc  length 
(Equation  ( 11-37. C)  are  found  directly.  When  constant 
elements  approximate  displacements,  finite  difference 
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methods  are  typically  used  to  provide  the  arclength 
derivatives  [65].  The  seven  equations  assembled  in  matrix 
form  appear  as 


1 

0 
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(11-38) 


where  the  traction  rates  and  the  elastic  tangential 
displacement  rate  derivatives  are  a product  of  the  boundary 
element  constraint  equation,  A = 2Gv/(l-2v),  a = -\-2v  , ni 
are  the  normals  and  G is  the  shear  modulus.  The  equations 
for  stress  rates  are  calculated  with  Cramer's  rule  and 
listed  below: 


an  = 


tlA  + G 4-  Gn|  ClC2 


D 


; (11-39. A) 


22  = 


-GB2E  + nj  Cl(t2n2-G  C2) 


(11-39. B) 


D 


tl  E + n2CiC2  + Vl 


ai2  = Gni( 


D 


} 


(11-39. C) 
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where 


A = A2n2+ni-a2nin22 , 

B1  = an12+An22,  B2  = n1t1-n2t2, 

• • 

Cl  = A2-a2,  C2  = u2,snru1;Sn2, 

D = Gan^+(2GA  -a2-A 2 ) ni2n2+Gan^ , 

E = Ani2+an2. 

The  final  point  of  numerical  interest  to  be  mentioned 
revolves  around  the  discontinuous  normals  at  corners.  The 
non-unique  normals  at  surface  discontinuities  lead  to 
discontinuous  surface  tractions.  The  implementational 
problems  incurred  by  linear  BEM  are  best  described  by 
example.  Referring  to  Figure  10,  both  element  A+  and 
element  A-  rightly  claim  node  A,  but  neither  "knows"  the 
other  has  staked  a claim.  The  problem  appears  when  elements 
A-  and  A+  both  try  to  assign  a value  to  the  tractions  at 
node  A based  upon  their  individual  normals  (n^+  and  n^_) . 

The  result  is  hard  to  predict  but  is  usually  erroneous  [56]. 
Other  than  physically  rounding  the  corners,  the 
discontinuity  problem  is  handled  by  either  the  double  node 
method  [17]  or  Chaudonneret ' s method  [66],  The  numerical 
hazards  of  lack  of  uniqueness  are  an  artifact  of  either 
inconsistent  linear  models  or  higher  order  models  in 
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general.  In  this  BEM,  the  tractions  are  constant  over  the 
element  and  located  at  the  element's  mid  point;  hence,  no 
elements  share  tractions.  Therefore,  each  element  has  a 
unique  traction  node  and  no  competition  occurs.  The 
discussion  of  the  surface  discontinuity  is  included  partly 
because  the  more  prevalent  inconsistent  BEM  models  encounter 
this  problem,  so  a discussion  of  the  non-unique  normals  is 
expected.  However,  the  comments  are  made  chiefly  to 
emphasize  the  advantage  afforded  by  this  consistent 
approach. 


NODE  A 


Figure  10.  Corner  node  with  a discontinuous 
normal . 

With  the  preliminary  numerical  eccentricities  taken 
care  of,  it  is  now  time  to  present  the  final  part  of  the 
iterative  EPBEM  solution  puzzle.  First  the  boundary 
integral  equation  for  the  total  strain  rate  is  needed.  By 
substituting  Equation  (11-27. A)  into  strain-displacement 
equation,  the  total  strain  rate  looks  generically  like 
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ui,k  = tj (s) 


Uij,k(s/P)dr  " 


Uj  (s)Tij /k(s'P)dr 


(11-40) 


+ 6 P-ii  d- 


D™ 


3xk 


^ijxn(s'P) 


where 

uij,k  = ij  (3-4^ ) r,  k “ 5jkr,i  + 2r,ir,jr,k) 

8ttG  ( 1—v  ) r 

7 (2  [ (2i/-l)  5 ijr  ^ + 5 j^r  i 
Jr  ( l—i/ ) r2 

r,ir,jr,k]  (— ) + (1—21/)  [nk5ij 

n 

i5jk  -2nir/jr/k  + 2njr/ir/k]  + 2nkrfirfj), 

* . • 
sijm  (s,p)dn  is  given  in  appendix. 

As  previously  mentioned,  the  singular  nature  of  Sjki(s,p) 
prevents  bringing  the  derivative  directly  under  the 
integral;  therefore  Sjki(s,p)  is  analytically  integrated  and 
then  differentiated  (Appendix) . 


Tij  /k  = 


- 4 


+ r 


and 


2_ 

axk 
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The  iterative  elasto-plastic  boundary  element  solution 
technique  is  similar  to  the  successive  approximation  scheme 
used  by  Roberts  and  Mendelson  for  stress  function  solutions 
[67].  The  iterative  approach  used  here  is  less  cumbersome 
since  the  total  strain  rates  are  found  directly.  The 
instructions  below  list  the  steps  in  the  EPBEM  solution 
technique  and  Figure  11  provides  a flow  chart  for  quick 
referencing. 


1)  Apply  an  elastic  load  to  bring  the  most 
highly  stressed  cell  to  a stress  state  near 
yield. 

2)  Find  the  total  elastic  stress  and  strain  at 
each  cell.  For  each  cell,  use  the  elastic 
strain  as  the  initial  guess  for  accumulated 

and  use  elastic  stress  as  the  initial 
guess  for  the  current  total  stress. 

3)  Apply  a load  increment  large  enough  to  cause 
incipient  yielding.  Use  the  load  increment 
and  the  initial  guess  for  the  plastic  strain 
rate  in  the  following  form  of  Equation  (II- 
31)  to  find  the  unknown  boundary  conditions. 

[A] { v}=  {q}  + 8 (q)  + [AS]{?p) 

5{q}  = load  increment 

4)  Use  {v}  from  step  (3)  to  find  the  total 
strain  rates  at  each  cell  with 
Equation  (11-40) . 

5)  Calculate  the  elastic  strain  rate  at  each 
cell  with  «eij  = «ij  ~ €pij. 

6)  Use  Hooke's  Law  (11-12)  to  calculate  the 
stress  rate  at  each  cell. 


Add  the  stress  increment  to  the  current  total 
stress  in  each  cell  and  use  it  to  find  the 
corresponding  tangent  modulus  (E-t)  from  the 
uniaxial  stress  strain  curve. 


V) 
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8)  Check  each  cell  to  see  if  it  is  beyond  yield; 

CTef  >?  Y'* 

where  a ef  is  the  von  Mises  equivalent  stress 
(also  the  yield  function)  and  Y is  the 
current  uniaxial  yield  strength  for  the  cell. 

9Y)  If  the  cell  has  yielded;  then  calculate  the 
new  plastic  strain  rate  of  the  cell  with  the 
isotropic  work  hardening,  von  Mises  flow  rule 
(11-33. C).  Use  the  new  total  stress,  Et  from 
step  (7)  and  the  current  plastic  strain  rate 
in  Equation  (11-33. C). 

9N)  If  the  cell  is  still  elastic,  set  ePfj  equal 
to  zero. 

10)  Check  for  ePfj  convergence  against  the 
previous  value. 

UN)  If  the  tolerance  is  not  met,  subtract  the 
stress  rate  found  in  step  (6)  from  the 
current  total  stress  found  in  step  (7) . Then 
go  to  step  (5)  and  start  again. 

11Y)  If  the  error  tolerance  is  met,  approximate 
the  new  uniaxial  yield  stress  for  each 
yielded  cell  by  the  current  von  Mises 
equivalent  stress. 

12)  If  the  problem  is  done,  quit.  If  the  problem 
is  not  finished  go  to  step  (3)  and  start  a 
new  increment  cycle. 


With  possible  exception  of  step  (11N) , the  instructions 
listed  above  are  self-explanatory.  For  step  (11N) , the 
current  guess  for  the  stress  rate  must  be  subtracted  from 
the  known  previous  total  stress  after  each  pass  through  an 
unsuccessful  iteration.  It  is  improper  to  add  each  current 
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Figure  11.  Flow  Chart  of  the  iterative  elasto-plastic 
boundary  element  solution  technique. 
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guess  for  the  stress  rate  to  the  last  increments  total 
stress  plus  the  last  iterations  stress  rate  guess.  The 
current  total  stress  would  then  be  the  sum  of  the  previous 
total  stress  and  all  the  stress  rate  guesses. 

The  above  iterative  solution  to  elasto-plastic  boundary 
element  problems  is  reasonably  efficient  and  has  acceptable 
accuracy.  The  technique  is  sufficiently  general  so  that  it 
is  easily  extended  to  other  than  boundary  element  methods. 
When  large  scale  yielding  occurs,  there  is  a possibility 
that  the  plastic  zone  will  exceed  the  celled  region.  If 
such  a case  arises,  one  expects  a reduction  in  EPBEM 
accuracy. 


II. 5 Numerical  Verification  of  EPBEM 

II. 5. a Plane  Strain  Circular  Cylinder  Under  Internal 

Pressure 


Two  well  documented  engineering  problems  in  plasticity 
serve  to  verify  and  evaluate  the  usefulness  of  the  initial 
strain  elastic-plastic  boundary  element  method  just 
described.  In  an  effort  to  isolate  the  initial  strain 
solution  approach,  the  plane  strain  radial  expansion  of  a 
long,  elastic-perfectly  plastic,  circular  cylinder  under 
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Figure  12.  Infinite  circular  cylinder  under 
internal  pressure. 
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internal  pressure  is  solved  by  EPBEM.  The  material 
behavior  is  approximated  with  the  von  Mises,  elastic- 
perfectly  plastic  model  and  the  numerical  results  are 
compared  to  the  solutions  presented  by  Hodge  and  White  [29]. 
The  purpose  of  this  example  is  to  show  the  validity  of  the 
successive  approximation  boundary  element  solution  scheme, 
independent  of  any  errors  introduced  by  the  work  hardening 
equations.  Figure  12  shows  the  geometry  and  loading 


Figure  13.  Domain  discretization  for  the  infinite 

circular  cylinder  under  internal  pressure. 
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conditions  for  the  cylinder.  Symmetry  enables  only  a 
quarter  of  the  cylinder  to  be  modelled.  Therefore,  the 
solution  requires  20  boundary  elements  and  25  internal  cells 
(Figure  13) . The  results  presented  in  Figure  14  and  Figure 
15  show  excellent  agreement  with  analytical  solutions  [29]. 
The  one  noticeable  deviation  is  located  at  the  elastic- 
plastic  boundary  (a=1.5r). 

The  small  deviation  is  assumed  to  be  a result  of  the 
algorithm  being  confused,  i.e.  not  knowing  whether  to 
consider  this  region  plastic  or  elastic.  Banerjee  et  al. 

[68]  evaluated  the  internal  stresses  near  the  elastic- 
plastic  interface  and  experienced  similar  deviation  in  the 
hoop  stress  but  not  the  radial  stress.  A discussion  of  this 
effect  is  not  presented  by  previous  authors  and  may  be 
restricted  to  the  present  algorithm.  More  complicated 
numerical  schemes  to  discriminate  a yielded  cell  did  not 
improve  the  results  enough  to  warrant  their  cost  in  CPU 
time. 

A major  advantage  of  elasto-plastic  boundary  element 
methods  lies  in  only  "zoning  up"  the  plastic  region. 

Placing  cells  only  in  the  plastic  region  saves  time  in 
setting  up  and  running  a problem.  This  advantage  can 
quickly  turn  into  a disadvantage  if  one  is  not  careful. 
Figures  16  and  17  show  the  results  of  the  circular  cylinder 
loaded  so  that  the  plastic  region  extends  to  1.6r.  In  this 
case,  the  solution  for  a <=  1.5r  is  correct  but  the 
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Figure  14.  Radial  Stress  Plot  for  the  infinite 
cylinder  under  internal  pressure. 


Figure  15.  Hoop  stress  plot  for  the  infinite 
circular  cylinder  under  internal 
pressure. 


NON-DIMENSIOAL  HOOP  STRESS  (St/Y)  NON-DIMENSIONAL  RADIAL  STRESS  (Sr/Y) 
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Figure  16.  Radial  stress  when  the  plastic  zone 
exceeds  the  celled  region. 


Figure  17 


Hoop  stress  solution  when  the  plastic  zone 
exceeds  the  celled  region. 
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solution  for  a =>  1.5r  goes  awry.  The  plastic  region 
exceeds  the  zoned  area  and  the  unaccounted  for  plastic 
deformation  results  in  appreciable  errors.  The  limited 
zoning  is  a great  advantage  but  it  is  also  perhaps  the  worst 
disadvantage  of  elasto-plastic  boundary  element  methods. 


Figure  18.  Infinite  plate  with  a central  hole. 
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II. 5. b Infinite  Uniaxial  Plate  With  a Central  Hole 

The  second  example  solved  by  EPBEM  includes  the  full 
working  hardening  equations  with  the  initial  strain  solution 
procedure.  An  infinite  plate  with  a central  hole,  loaded  as 
shown  in  Figure  18,  provides  sufficient  complexity  to 
challenge  the  EPBEM  routine.  This  geometry  is  modelled  in 
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Figure  19.  Model  region  of  the  infinite  plate  with  a 
central  hole. 
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Y 


Figure  20.  Domain  discretization  for  the  infinite 
plate  with  a central  hole. 


Bi-linear  uniaxial  stress-strain  curve 
used  in  the  solution  of  the  infinite 
plate  with  a central  hole. 


Figure  21. 
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much  the  same  way  as  Riccardella  [17].  The  model  of  a 
circular  region  of  10R  enclosing  the  cutout  is  sketched  in 
Figure  19.  In  addition,  the  modelled  plastic  zone  is  shown 
in  Figure  20  and  material  properties  are  given  in  Figure  21. 

Two  cases  are  explored.  Both  cases  are  compared  to  the 
elastic  solution  loaded  so  that  the  most  highly  stressed 
region  is  at  the  initial  yield  point  [50].  They  are  also 
compared  to  Swedlow's  isotropic  hardening,  finite  element 
solution  [69]  (as  reported  by  Riccardella  [17]).  In  the 
first  case,  loading  continues  until  the  ratio  of  the  applied 
stress  to  the  current  yield  stress  is  0.4;  this  ratio  is 
known  as  the  load  factor  (L.F.).  Incipient  yielding  occurs 
at  a load  factor  of  0.333  [50];  therefore  the  L.F.=0.4 
results  should  not  deviate  appreciably  from  the  elastic 
results.  This  observation  is  demonstrated  by  Figure  22. 
Figure  23  shows  the  equivalent  stress  concentration  factors 
for  a load  factor  of  0.8.  As  expected  the  results  do  not 
correlate  with  the  elastic  solution  but  do  show  good 
agreement  with  Swedlow's  results. 


STRESS  CONCENTRATION  FACTOR  (St/So)  STRESS  CONCENTRATION  FACTOR  (St/So) 
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NON-DIMENSIONAL  DISTANCE  (X/R) 

ELASTIC  X FEM  V BEM 

Figure  22.  Stress  concentration  factor  for  a 
load  factor  of  0.4 


Figure  23.  Stress  concentration  factor  for  a 
load  factor  of  0.8. 
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II. 6 Parametric  Study  of  Experimental 
Effects  on  the  EPBEM  Algorithm 

This  final  section  explores  the  influence  experimental 
errors  have  on  the  accuracy  of  a hybrid  experimental- 
elastic-plastic  boundary  element  method.  Sutton  et.  al . 
[70]  emphasizes  that  errors  in  measuring  boundary  data  can 
sometimes  have  a very  large  impact  on  solution  on  or  near 
the  boundary.  Heeding  these  comments,  the  following  three 
types  of  experimental  error  will  be  studied  from  a 
parametric  viewpoint:  1)  system  error,  2)  large  error  at 

one  data  point,  and  3)  both  system  error  and  large  error  at 
one  data  point.  To  demonstrate  clearly  the  induced  errors, 
each  test  explores  the  response  of  a unit  square  under 
uniaxial  tension  (Figure  24)  with  the  material  properties 
given  in  Figure  21.  For  completeness,  the  y-traction  at 
mid-node  2 and  the  y-strain  at  cell  4 are  monitored  and 
compared  to  results  of  a similar  completely  elastic  unit 
square. 

Every  cell  in  the  unit  square  yields  at  the  same  time. 
This  eliminates  any  errors  introduced  by  the  elastic-plastic 
interface  and  also  eliminates  any  redeeming  features 
resulting  from  local  yielding.  If  a large  plate  with  a 
small  hole  under  the  same  loading  conditions  were  used  to 
evaluate  the  influence  of  experimental  boundary  data  error 
on  the  interior  stress  state,  one  would  not  expect 
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Figure  24.  Unit  square  used  in  the  parametric  study. 

perturbations  in  boundary  data  to  greatly  influence  the 
plastic  stress  state  in  the  region  localized  near  the  hole. 
This  is  a physical  argument  based  on  geometry  and  the 
loading  conditions.  The  purpose  of  this  parametric  study  is 
to  evaluate  the  numerical  stability  of  the  internal  stresses 
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and  strains  in  relation  to  boundary  perturbations. 

Therefore,  the  unit  square  is  the  ideal  choice  to  eliminate 
local  yielding  situations. 

The  first  test  will  explore  system  error  effects  on  the 
EPBEM  solution.  System  error  refers  to  general  errors 
inherent  to  the  equipment,  experimental  technique, 
environmental  considerations,  etc.  Random  and  zero  mean 
errors,  ranging  from  0%  to  20%  are  added  to  the  exact 
boundary  conditions  shown  in  Figure  24.  The  results  of 
Figure  25  reveal  that  the  internal  plastic  strain  solution 
is  quite  accurate  until  the  random  error  reaches  16%.  Then, 
the  solution  appears  to  diverge  from  the  correct  result. 

Even  so,  the  plastic  strain  solution  is  much  less  sensitive 
to  boundary  perturbation  than  the  elastic  strain  solution 
for  an  equivalent  elastic  unit  square.  On  the  other  hand, 
the  elastic  and  plastic  responses  at  mid-node  2 are  very 
similar  (Figure  26) , both  diverging  immediately. 

When  an  experimentalist  finds  a bad  data  point  he 
usually  just  disregards  it  and  tries  again.  Automated  data 
extraction  techniques  may  not  have  this  ability. 

Consequently  the  influence  of  a single  bad  data  point  on  the 
accuracy  of  the  EPBEM  must  be  explored.  In  this  example,  a 
range  of  0%  to  20%  error  is  added  to  the  y-traction  at  mid- 
node 7 (Figures  27  and  29)  and  then  20xl0-6m  (in  lxl0-6m 
increments)  is  added  to  the  x-displacement  at  node  4 
(Figures  28  and  30).  In  these  figures,  the  internal  plastic 
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strains  always  remain  well  within  acceptable  error  for  large 
perturbations  but  the  tractions  immediately  starts  to 
diverge. 

A series  of  plots  showing  the  influence  of  both  system 
and  single  point  error  is  given  in  Figures  31  through  34. 
These  plots  show  only  plastic  response  since  the  elastic 
response  is  linear  and  are  inferible  from  the  previous 
graphs.  The  nonlinearity  of  the  problem  is  evident  in 
Figures  31  and  33,  where  the  internal  strain  starts  to 
diverge  at  a much  earlier  error  in  the  y-traction  at  mid 
node  7 than  might  be  predicted  for  the  earlier  graphs. 

What  can  be  inferred  from  all  this  error  analysis? 
Except  for  a feel  for  acceptable  error  limits,  little 
guantitative  can  be  gained.  For  instance,  the  plastic 
strain  is  much  less  sensitive  to  boundary  data  error  than 
its  elastic  counterpart.  Coupled  with  local  yielding 
arguments,  this  information  can  be  comforting.  On  the  other 
hand,  the  elastic-plastic  and  purely  elastic  tractions  are 
egually  sensitive  to  surface  condition  perturbations.  Since 
the  most  interesting  region  is  normally  located  near  a 
boundary,  the  parametric  error  study  tells  the  researcher 
what  he  already  knows:  be  meticulous  when  searching  for 
information  near  a boundary.  Excellent  accuracy  in  high 
stress  concentration  regions  is  attainable  [52,  71-72]  but 
one  must  be  aware  of  the  large  influence  of  experimental 
error  in  order  to  properly  evaluate  EPBEM  results. 


Y-TRACTION  AT  MID-NODE  2 PERCENT  ERROR  IN  Y-STRAIN  AT  CELL 
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PERCENT  RANDOM  ERROR  IN  B.C. 
n ELASTIC  BEM  V PLASTIC  BEM 


Figure  25.  Internal  elastic  and  plastic  strain 
response  to  random,  zero  mean  error 
added  to  the  exact  boundary  conditions. 


PERCENT  RANDOM  ERROR  IN  B.  C. 
n ELASTIC  BEM  V PLASTIC  BEM 

Figure  26.  Elastic  and  plastic  traction  response  to 

random,  zero  mean  error  added  to  the  exact 
boundary  conditions. 


Y-TRACTION  AT  MID-NODE  2 PERCENT  ERROR  SN  Y— STRAIN  AT  CELL  4 
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% ERROR  IN  Y-TRACTION  AT  MID-NODE  2 
n ELASTIC  BEKI  V PLASTIC  BEM 


Figure  27.  Elastic  and  plastic  strain  response  to 
errors  in  traction  at  one  location. 


X ERROR  IN  Y-TRACT.ON  AT  MID-NODE  7 
D ELASTIC  BEM  V PLASTIC  BEM 

Figure  28.  Elastic  and  plastic  traction  response  to 
to  error  in  traction  at  one  location. 


Y-TRACTLON  AT  MID-NODE  2 PRECENT  ERROR  IN  Y— STRAIN  AT  CELL  4 
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X— DISPLACEMENT  AT  NODE  4 (1 0X-6m) 
n BASTiC  BEM  V PLASTIC  BEM 


Figure  29.  Elastic  and  plastic  strain  response  to 

error  added  to  displacement  at  one  point. 


X— DISPLACEMENT  AT  NODE  4 (10X-6m) 
n ELASTIC  BEM  V PLASTIC  BEM 


Figure  30.  Elastic  and  plastic  traction  response  to 
error  added  to  displacement  at  one 
location. 


•TRACTION  AT  MID-NODE  2 PERCENT  ERROR  IN  Y-STRAIN  AT  CELL  4 
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SJC  RANDOM  ERROR  IN  B.  C. 


Figure  31.  Plastic  strain  response  to  both 
random  error  and  single  location 
traction  error. 


5Z  RANDOM  ERROR  IN  B.  C. 


Figure  32.  Plastic  traction  response  to  both  random 
error  and  single  location  traction  error. 


Y— TRACTION  AT  MID-NODE  2 % ERROR  IN  STRAIN  AT  MID-NODE  4 
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5%  RANDOM  ERROR  IN  B.  C. 


Figure  33.  Plastic  strain  response  to  both  random 
error  and  single  location  displacement 
error. 


5 % RANDOM  ERROR  IN  B.  C. 


Figure  34.  Plastic  traction  response  to  both  random 
error  and  single  location  displacement 
error. 


CHAPTER  III 


DISPLACEMENT  PATTERN  MATCHING 


The  spread  of  image  processing  applications  in 
photomechanics  is  a consequence  of  the  availability  of 
equipment  at  affordable  costs.  Hardware  prices  continue  to 
decline  and  the  power  of  personal  computers  continues  to 
increase.  Economical  mass  storage,  networking,  graphics, 
and  application  specific  software  all  have  increased  the 
proliferation  of  image  processing  in  the  photomechanics 
field  [73].  For  the  most  part,  image  processing  techniques 
in  photomechanics  are  grouped  into  two  categories,  fringe 
analysis  [74-77]  and  correlation  techniques  [6,78-80].  A 
third  category,  pattern  mapping  [23,81],  is  emerging  as 
alternative  with  many  attractive  qualities.  Pattern  mapping 
is  a non-destructive  technique  measuring  displacement  and 
strain  and  is  based  on  image  processing  and  syntactic 
pattern  recognition.  It  is  a process  by  which  "black"  spots 
on  a "white"  background  are  followed  in  a double  exposure 
type  scheme.  This  technique  offers 
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1)  Full  field  measurement. 

2)  Sub-pixel  registration. 

3)  Greatly  reduced  CPU  time. 

4)  A highly  automated  technigue. 

5)  Relative  ease  of  specimen  preparation. 

Pattern  mapping  was  originally  developed  as  a general 
purpose  strain  measurement  technique  capable  of  discerning 
large  rigid  body  rotation  and  translations.  It  also 
included  a variety  of  strain  definitions.  Fail  [23] 
suggested  that  reductions  in  run  time  were  achievable  by 
tailoring  the  pattern  mapping  to  a specific  class  of 
problems.  Displacement  pattern  matching  ( DPM)  is  a result 
of  this  suggestion.  The  displacement  pattern  matching 
algorithm  is  developed  by  assuming  that  the  experimenter  can 
prevent  appreciable  rigid  body  motion  and  that  displacements 
are  the  only  information  of  interest.  As  a result,  a highly 
automated,  very  fast  displacement  tracking  scheme  emerged. 

The  DPM  algorithm  can  be  broken  into  three  blocks: 
image  segmentation,  feature  extraction,  and  matching. 

These  three  terms  are  common  pattern  recognition  jargon  but 
may  be  new  to  the  experimental  mechanics  community.  To 
clarify  the  transition  into  the  new  terminology,  formal 
definitions  of  image  segmentation,  feature  extraction, 
matching  and  a fourth  term,  registration,  are  provided. 
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Image  Segmentation  is  the  separation  of  the  image 
into  different  regions,  each  having  specific  and 
related  properties.  It  is  the  first  step  in  all 
pattern  recognition  schemes  and  attempts  to 
describe  or  classify  the  image  [82,83].  An 
example  might  be  to  describe  a landscape;  that  is, 
segmentation  of  a satellite  photograph,  where 
classification  might  entail  distinguishing  rivers 
from  the  background  [84]. 

Feature  Extraction  is  the  quantification  of 
certain  characteristics  of  an  image  which  are 
usable  in  further  image  information  classification 
[85].  This  process  discerns  important  foreground 
features  from  the  background.  Examples  might  be 
lengths,  areas,  edges  or  curves  [86]. 

Matching  refers  to  a class  of  operations  of 
comparing  to  two  or  more  images  [87].  Matching  is 
typically  used  in  time-varying  imagery,  where 
motion  detection  is  simply  a matter  of  detecting 
differences  in  successive  images  [88,89]. 

Registration  is  the  comparison  of  two  images  of  a 
scene  taken  from  different  perspectives. 
Registration  is  used  as  a measure  of  resolution 
attainable  with  matching  algorithms  [90] . 


Capturing  a digital  image  and  the  equipment  required  to  do 
so  are  briefly  discussed  in  the  next  section.  The  ensuing 
sections  describe  the  segmentation,  feature  extraction,  and 
matching  blocks  of  the  DPM  algorithm.  The  final  section 
verifies  DPM  through  rigid  body  motion  tests  and  will 
attempt  to  assess  the  important  characteristics  that  control 


accuracy. 
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III.l  Images  and  Digitizing  Equipment 

The  human  eye  scans  the  analog  intensity  distribution 
which  is  created  by  light  reflecting  from  a scene  and  sends 
this  information  to  the  brain  to  be  interpreted.  Upon 
receiving  the  signal,  the  brain  processes  the  scene  and 
classifies  the  image.  Ideally,  computer  vision  should  be 
capable  of  emulating  this  process.  A quantization  of  the 
intensity  distribution  is  necessary  if  a computer  is  to  have 
any  hope  of  mimicking  the  recognition  processes  of  the  human 
brain.  The  most  primitive  medium  of  communication  for  a 
computer  is  numbers  (digits) . Therefore,  it  is  logical  that 
an  efficient  quantization  technique  would  be  able  to  convert 
the  analog  intensity  distribution  to  a digital 
approximation,  hence  the  term  digitize. 

A digital  image  is  a numerical  approximation  of  a real 
scene,  with  local  intensity  distributions  approximated  by 
gray  levels  [91] . Gray  levels  are  integers  ranging  from 
zero  to  an  upper  limit  set  by  the  hardware.  The  eight  bit 
digitizing  board  used  in  these  experiments  affords  0-255 
gray  levels.  These  gray  levels  represent  the  darkness  and 
are  ordered  to  form  a matrix  representation  of  the  image. 
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An  example  of  the  imaging  procedure  is  given  in  Figure  35 
where  a black  spot  and  its  digital  representation  are 
provided.  Key  traits  to  note  include:  black  appears  as 
lower  gray  levels  than  white,  the  two  pixel  (picture 
element)  fuzzy  transition  region  between  black  and  white 
[92],  and  the  noise  in  the  image.  When  looking  at  scenes 
containing  high  contrast  sub-regions,  the  human  eye  has 
little  difficulty  locating  boundaries  between  light  and 
dark.  For  computer  vision  hardware,  this  process  is  not  so 
simple.  One  would  expect  sharp  contrast  borders  in  the 
scene  to  be  represented  by  a jump  in  gray  level,  when  in 
fact,  the  imaging  system  sees  a finite  transition  between 
light  and  dark.  This  fuzzy  transition  is  a product  of 
hardware  limitations  in  the  camera  and  in  the  digitizing 
board.  The  camera  sensing  array  is  divided  into  a finite 
number  of  sensing  cells.  Cell  crosstalk,  frequency 
limitation  characteristics,  and  cell  overlap  are  typical 
contributions  of  the  camera  to  fuzzy  regions  [93-94].  in 
addition,  the  mismatch  between  the  camera  array  size  and  the 
digitizing  hardware  array  size  increases  the  fuzziness.  In 
general  these  array  sizes  differ  considerably.  A typical 
example  would  be  four  digitizing  pixels  representing  three 
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Figure  35.  A black  spot  and  its  digital 
representation . 
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camera  pixels,  i.e.  a smearing  effect.  Although  the 
mismatch  is  not  usually  this  dramatic,  its  effect  is 
noticeable. 

In  Figure  35  the  spot  is  uniformly  black;  unfortunately 
noise  degradation  prevents  the  image  from  reflecting  the 
uniformity.  Degradation  is  apparent  in  the  fluctuating  gray 
level  representation  of  black  in  the  spot  and  results  from 
optical  or  electronic  sources.  Spherical  aberration,  coma, 
astigmatism,  non-uniform  illumination,  and  dust  are  just  a 
few  examples  of  optical  noise  [95].  Electronic  noise  may 
come  from  thermal  effects  in  the  imaging  array,  the 
digitizing  process,  electro-magnetic  interference,  or  may 
just  be  inherent  to  the  imaging  array  [92-94].  Regardless 
of  the  source,  DPM  regards  noise  as  random  and  therefore 
having  a zero  mean  distribution. 

Some  comments  about  the  image  sampling  are  now 
appropriate.  The  process  of  digitizing  an  object  involves 
sampling  and  quantizing  the  intensity  distribution  at  IxJ 
points  (sampling  spatially  and  in  amplitude) . The  exact 
representation  of  a continuous  field  requires  an  infinite 
number  of  pixels  in  IxJ  image  matrix  and  an  infinite  number 
of  gray  levels.  Since  infinite  limits  are  impossible, 
limits  must  be  established.  This  section  touches  on  several 
important  consequences  of  limiting  the  number  of  gray  levels 
and  sample  points  [96-98]. 
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Figure  36.  An  image  of  a lady's  face  sampled 

with  384H  x 485V  spatial  resolution 
and  0-255  gray  scale. 
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Figure  37 


An  image  of  the  lady's  face  with 
96H  x 120V  spatial  resolution  and 
0-255  gray  scale. 
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Digital  image  resolution  is  highly  dependent 
on  the  number  of  spatial  sample  points.  When  considering 
the  extreme  example  of  representing  an  entire  scene  with  one 
sample  point,  it  is  apparent  that  the  result  in  no  way 
resembles  the  scene.  This  suggests  a lower  bound  on  the 
number  of  sample  points  reguired  to  represent  an  object 
adeguately.  Undersampling  of  an  image  has  the  same  effect 
as  a low  pass  filter;  the  image  can  become  fuzzy  to  the 
point  that  it  is  not  recognizable  [92].  Figures  36  and  37 
show  a digital  image  of  a lady  portrayed  with  both  384H  x 
485V  and  96H  X 120H  pixels.  Note  the  fuzzy  lines  of 
contrast  in  the  face,  the  greatly  reduced  resolution  of  her 
ear,  and  the  aliasing  of  the  Venetian  blinds  in  the  96H  x 
120V  representation.  Image  processing  requires  good 
resolution  to  be  successful;  Gonzalez  and  Wintz  [89]  suggest 
a matrix  size  of  256H  x 256V  as  the  minimum  size  needed  to 
produce  an  image  of  adequate  quality. 

In  contrast  to  undersampling,  oversampling  is 
defined  more  by  storage  space  and  field  of  view  than  by 
image  representation.  Ideally,  one  would  like  to  use  as 
many  sample  points  as  conceivable,  but  practically,  image 
storage  size  severely  taxes  the  computer  storage  capacity. 

A 512H  x 512V  one  byte  element  array  requires  more  than  a 
quarter  of  a megabyte  of  storage.  Further,  image  processing 
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Figure  38.  An  image  of  the  lady's  face  with 

384H  x 485V  spatial  resolution  and  0-3  gray 
scale. 
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and  pattern  recognition  are,  on  the  whole,  computationally 
expensive,  run  time  increases  dramatically  with  the  sample 
set  size. 

If  one  desires  very  high  resolution,  a zoom  lens  can 
localize  the  imaging  array  on  a small  portion  of  the  scene. 
In  this  case,  many  pixels  digitize  one  part  of  a scene  but 
the  expense  is  not  sampling  other  parts  of  the  scene  at  all. 
The  cost  is  a severely  restricted  field  of  view  resulting  in 
a significant  loss  of  information  [89,97]. 

Finally,  amplitude  quantization  needs  to  be  explored. 
Gray  scale  is  the  range  of  possible  gray  levels  available  to 
represent  a scene  and  is  hardware  defined.  A simple 
experiment  can  illustrate  the  effect  of  gray  scale  on  a 
digital  image.  The  lady's  face  in  Figure  36  (0-255  gray 
levels)  is  now  digitized  using  a gray  scale  0-3  (2  bits)  and 
shown  in  Figure  38.  False  contouring  and  poor  resolution 
degrade  the  image.  Degradation  is  particularly  noticeable 
in  the  poor  background  resolution. 

Now  that  the  ideas  and  important  aspects  of  a digital 
image  are  illustrated,  equipment  and  procedure  for  obtaining 
a digital  image  can  be  described.  Figure  39  depicts  a 
typical  personal  computer  based  digitizing  system.  The 
imaging  device  used  here  is  a Sony  XC-38  CCD  camera,  a 
Datacube  IVG128  Video  Acquisition  and  Display  Board 
(384H  x 512V)  which  quantizes  the  image  with  its  digital 
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VIDEO  MONITOR 


CCD  CAMERA 


Figure  39.  Typical  personal  computer  based  digitizing 
system. 


form  displayed  on  a PVM-1271Q/1371QM  Sony  studio  monitor. 
All  image  processing  is  done  on  a PC's  Limited  286-8 
Personal  computer  equipped  with  a math  co-processor  and  EGA 
graphics . 

Digitizing  is  software  controlled.  On  command,  the 
frame  buffer  stores  the  image  as  a matrix  of  signed  bytes 
representing  the  gray  levels.  This  two-dimensional  matrix 
is  also  available  to  the  user  on  software  control.  On 
software  control,  the  analog  signal  from  the  CCD  (charge 
couple  device)  camera  is  quantized  with  an  analog-to-digital 
converter  and  then  mapped  to  the  appropriated  gray  levels 
by  a software  selectable  look-up  table.  The  image  size 
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stored  in  the  frame  buffer  is  384H  x 485V  to  conform  with 
the  RS-170A,  60Hz  United  States  interlaced  video  standard. 
The  3 84H  x 512V  image  conforms  to  the  CCIR,  50Hz.  European 
interlaced  video  standard.  The  frame  buffer  is  accessible 
to  the  CPU  (central  processing  unit)  through  the  personal 
computers  I/O  page  [99].  VNA  Systems  Incorporated's  Digital 
Correlation  Metrology  software  package  provides  software 
control  of  the  digitizing  process. 


III.  2 Segmentation 

A decided  advantage  of  displacement  pattern  matching  is 
the  limited  prior  knowledge  required  for  success.  Success 
of  the  technique  hinges  on  the  placement  of  high  contrast 
spots  (dark  spots  on  a light  background)  at  locations  where 
displacement  is  desired.  Implicit  in  this  knowledge  is  the 
idea  that  the  images  are  binary,  that  is,  the  gray  levels 
should  be  either  black  or  white  [100].  In  keeping  with  the 
automated  methodology,  DPM  does  not  need  prior  knowledge 
about  the  number  of  spots,  spot  area  or  any  other 
distinguishing  feature.  The  primary  function  of 
segmentation  is  to  locate  the  important  regions  of  an  image. 
Since  DPM  images  are  binary,  segmentation  consists  of 
identifying  those  individual  regions  defined  as  black  (the 
spots)  since  they  contain  the  desired  information  [101]. 


81 


Segmentation  typically  (as  is  done  here)  starts  with 
some  form  of  characteristic  feature  thresholding  [102].  In 
essence,  a value  of  some  characteristic  feature  of  the  image 
is  chosen  as  a basis  on  which  to  make  decisions.  This 
pivotal  value  is  known  as  the  "threshold".  In  DPM,  the 
characteristic  feature  is  the  gray  level  and  the  threshold 
is  the  parameter  used  to  discriminate  between  black  and 
white.  All  gray  levels  above  the  threshold  are  white  and 
all  gray  levels  below  are  black.  The  transformation  to  a 
true  binary  gray  scale  enables  a border  following  scheme  to 
locate  the  spots.  Discerning  if  each  spot  is  a true  spot  or 
if  it  should  be  discarded  is  segmentation's  culminating 
decision. 

Since  it  is  the  first  step  in  finding  spots,  threshold 
choosing  is  perhaps  the  most  critical  phase  of  DPM.  By 
definition,  general  thresholding  techniques  [102-104]  do  not 
use  information  inherent  to  the  image  to  approximate  the 
threshold.  These  techniques  are  at  most  over-deterministic 
and,  at  the  least,  computationally  expensive.  Estimating 
the  gray  level  threshold  through  basic  image  statistics 
offer  an  efficient  alternative  [105-106]. 

The  gray  level  frequency  distribution  or  "histogram," 
of  binary  images  will  always  be  bi-modal.  Only  two  spikes 
are  present  in  ideal  binary  images  (black  and  white) 


NUMBER  OF  PIXELS  AT  EACH  GRAY  LEVEL 
(Thousands) 
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but  the  aforementioned  noise  degrades  the  image.  Figure  40 
is  a histogram  of  the  binary  image  containing  the  array  of 
spots  shown  in  Figure  41.  As  seen  in  Figure  40,  the 
histogram  shows  two  groups  (or  modes) . One  is  centered  on 
black  and  the  other  is  centered  on  white;  hence  the  term  bi- 
modal.  By  noting  the  threshold  lies  somewhere  in  the  valley 
of  the  two  modes  [107],  a guick  estimation  of  the  gray  level 
threshold  is  possible. 


Figure  40.  Histogram  of  Figure  42. 
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Figure  41.  An  image  of  a 5x5  square  array  of 
spots. 
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The  threshold  estimation  process  begins  by  randomly 
sampling  sixty  pixels  and  then  calculates  the  average  gray 
level,  1^,  and  the  standard  deviation,  o&.  Next,  a mid- 
valley gray  level  value  is  defined  as  1^  = Ia-ctA  (see 
Figure  40) . Sixty  pixels  with  gray  levels  below  IM  are 
randomly  sampled  and  averaged  to  find  1 3 (average  black) . 
Sixty  more  pixels  with  gray  levels  above  1^  are  averaged  to 
find  Iw  (average  white) . The  threshold  (IT)  is  defined  by 
the  average  of  1^  and  I3; 

XM  = (XW  + Ib)/2  (HI-1) 

A more  detailed  discussion  of  this  technique  can  be  found 
in  Ref.  [23].  As  a point  of  interest,  other  thresholding 
schemes  specifically  designed  for  bi-modal  histograms  have 
been  developed  (see  Rosenfeld  and  Kak  [92]  pg.  62)  but  are 
computationally  intensive. 

The  next  important  segmentation  step  is  a process 
called  raster  scanning.  A raster  scan  is  nothing  more  than 
a row  by  row  or  column  by  column  evaluation  of  each  pixel 
for  some  prescribed  property;  in  this  case  to  find  the  first 
spot  (gray  level  < I-p  ) . When  DPM  locates  a black  pixel, 
the  gray  levels  of  its  eight  nearest  neighbors  are  checked. 
If  at  least  one  neighbor  is  black  then  DPM  has  located  a 
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spot.  If  none  of  its  neighbors  are  black,  then  DPM 
considers  the  center  pixel  as  a false  spot.  Once  a true 
spot  is  located,  a general  border  following  scheme  developed 
by  Rosenfeld  [108]  and  specialized  to  only  find  outside 
edges  [109,23]  'frames'  the  spot.  Framing  the  spot  is  the 
process  by  which  a complete  clockwise  circuit  around  the 
border  is  traversed,  storing  only  the  left  most,  right  most, 
highest,  and  lowest  pixel  locations  (L,  R,  T,  B 
respectively)  for  the  feature  extraction.  To  account  for 
the  transition  between  black  and  white,  the  frame 
automatically  expands  by  two  pixels  in  all  directions 
(Figure  35) . With  the  first  spot  framed,  the  raster 
scanning  continues.  This  process  is  repeated  until  the 
entire  image  is  examined. 

False  spots  larger  than  one  pixel  in  diameter  are 
entirely  possible,  especially  when  the  scene  is 
significantly  magnified.  An  automated  method  of  discerning 
false  from  true  spots  is  therefore  required.  Experience 
indicates  that  the  areas  of  false  spots  are  much  smaller 
than  their  legitimate  cousins  (area  is  defined  as 
A = [R-L][B-T]).  A spot  area  larger  than  one  third  the 
average  spot  area  (averaged  over  every  framed  spot)  is  a 
true  spot.  This  definition  is  completely  arbitrary  but  it 
proved  to  be  99%  reliable  for  the  spot  sizes  used  in  this 


research. 
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III. 3 Feature  Extraction 


Features  are  selected  measurements  which  are  either 
invariant  or  less  sensitive  to  typical  system  distortions 
[110];  they  tend  to  distinguish  the  object  from  its 
background.  Two  types  of  features  are  present  in  the  DPM 
images,  statistical  and  structural  [110].  The  statistical 
features  are  the  frame  boundaries:  L,  R,  T,  B,  and  the  spot 
centroid.  The  structural  features  are  the  spots  themselves. 
Since  segmentation  simultaneously  extracts  the  frame  and  the 
spots,  only  the  centroids  need  be  extracted. 

The  spot  centroids  are  calculated  using  the  true  gray 
levels  in  the  frame  and  pixel  locations  (i,j)  in  the  matrix 
[92].  Mathematically  stated  the  centroids  are 
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where  I(i,j)  is  the  gray  level  at  element  (i,j)  in  the  image 
matrix.  The  influence  of  the  zero  mean  noise  on  the 
centroid  calculations  is  presumed  to  average  to  zero. 


III. 4 Image  Matching 

Zero  rigid  body  motion  reduces  the  syntactic  pattern 
recognition  language  [111]  of  DPM  from  context  sensitive  to 
context  free  [112].  Pattern  mapping  reguires  each  spot  to 
be  located  at  the  same  "address"  before  and  after 
deformation  [23].  In  computational  space,  this  means  each 
spot  must  always  be  surrounded  by  the  same  neighbors.  If 
this  is  not  the  case,  that  spot  is  out  of  context  which 
means  it  makes  no  syntactic  sense.  DPM  makes  no  such 
requirement.  The  spot  locations  are  stored  in  the  order  in 
which  they  are  found.  Pattern  matching  [113]  is  then  used 
to  find  the  spots  in  all  ensuing  registered  images. 

Ordering  the  spots  by  their  correct  "addresses"  is  no  longer 
crucial  to  finding  the  next  image's  spots;  therefore,  the 
spots  do  not  need  to  be  in  context  to  make  syntactical 
sense. 

An  instructive  approach  to  describing  pattern  matching 
is  to  take  a figurative  step  back  and  simply  look  at  the 
physical  nature  of  a typical  solid  mechanics  problem.  The 
question  is  "what  are  the  displacements  at  points  a,  b,  and 
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c?"  To  answer  this  question  correctly,  each  spot  must  be 
examined  after  deformation  and  compared  to  its  previous 
location.  Successive  load  increments  produce  a sequence  of 
comparisons.  In  elasticity  and  incremental  plasticity, 
assuming  spot  motion  to  be  small  relative  to  spot  size  is 
restrictive  but  certainly  reasonable.  This  indicates  that 
each  sequential  image  should  closely  resemble  its 
predecessor.  The  above  description  of  events  is  not  unique; 
Rosenfeld  and  Kak  [108]  (pages  50-51)  give  an  almost 
identical  description  of  events  but  their  discussion 
pertains  to  time-varying  imagery. 

Time-varying  imagery  is  a popular  form  of  pattern 
matching  [87,114-115].  It  is  a sequential  registration 
technique  used  in  motion  detection,  object  tracking,  and 
dynamic  scene  analysis;  therefore,  it  is  well  suited  for 
discerning  displacements.  Here,  spot  motion  is  measured  by 
local  matching  of  the  segmented  spots.  In  essence,  the 
location  of  a spot  in  the  image  is  assumed  near  its  location 
in  the  previous  sequential  image.  Therefore,  instead  of 
raster  scanning  the  entire  new  image,  only  the  previous 
framed  area  is  scanned. 

If  the  displacement  field  carries  the  spot  out  of  the 
original  frame  area  then  the  search  region  is  expanded  to 


An  = [ L-R] [B-T] n2 , n=2 ,3,4...  (III-4) 
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n = 3 


Figure  42.  Search  region  expansion. 

where  An  is  the  area  of  the  search  region  and  n is  the 
number  of  search  region  expansions  required  to  find  the  spot 
(Figure  42) . When  the  search  region  encompasses 
more  than  one  spot  (Figure  43)  or  a wrong  spot  only 
(Figure  44) , a one-to-one  match  is  not  guaranteed  and  the 
matching  process  fails.  Therefore,  over-expansion  defines  a 
restriction  on  DPM,  but  at  the  same  time  somewhat  relaxes 
initial  assumptions  of  the  matching  scheme.  The  spot 
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SPOT  IN  THE 
PREVIOUS  IMAGE 


SPOT  IN  THE 
CURRENT  IMAGE 


n = 4 


Figure  43.  Search  region  encompassing  more  than  one 
spot. 


SPOT  IN  THE 
PREVIOUS  IMAGE 


SPOT  IN  THE 
CURRENT  IMAGE 


Figure  44.  Search  region  encompassing  the  wrong  spot. 
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displacement  must  be  smaller  than  or  equal  the  distance  from 
the  original  spot  to  the  closest  spot  in  the  current  image 
(Figure  44) . This  critical  distance  is  easily  gaged.  For 
example,  in  a square  array  of  spots,  the  maximum 
displacement  can  never  be  greater  than  0.7071  times  the 
smallest  original  spot  spacing.  Spot  displacements  are  no 
longer  restricted  to  being  small  compared  with  the  spot 
diameter. 

A special  search  region  expansion  case  occurs  near  the 
image  edge.  When  the  displacement  field  carries  a spot 
close  to  the  perimeter  of  the  image,  the  expansion  may  try 
to  extend  the  search  region  beyond  the  image  border,  the 
results  of  which  are  unpredictable.  The  search  region  is 
therefore  contoured  to  the  image  edge  if  the  region 
encounters  the  perimeter. 

With  a one-to-one  correspondence  established  between 
segments  in  successive  images,  statistical  feature 
extraction  can  proceed.  The  centroid  of  each  spot  is 
calculated  by  Equations  (III-2)  and  (III-3)  and  the  x-  and 
y-coordinates  of  the  centroids  in  sequential  images  are 
subtracted  to  find  the  displacements 

ux  = Xn-Xn_!  ( HI-5A) 

and 

uy  = Yn-Yn-1* 


(III-5B) 
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In  these  displacement  definitions,  n refers  to  the  current 
image.  The  units  of  displacement  are  pixels  but  can  be 
converted  to  a units  of  length  through  calibration. 

III. 5 Rigid  Body  Motion  Tests 

In  this  closing  section,  the  limitations  of 
displacement  pattern  matching  are  explored.  A series  of 
rigid  body  motion  tests,  each  with  different  spot  diameters, 
evaluate  the  accuracy  of  DPM.  In  the  process,  general 
accuracy  characteristics,  accuracy  versus  spot  size,  effect 
of  image  averaging  on  accuracy  and  the  consequences  of 
search  region  over-expansion  are  examined. 

Each  rigid  body  motion  test  consists  of  tracking  the 
5x5  array  of  spots  shown  in  Figure  42  through  four 
horizontal  0.025  (±0.001)  inch  increments.  Spot  areas  are 
easily  altered  by  magnification  with  a zoom  lens. 

These  tests  produce  results  that  are  difficult  to 
assimilate  and  make  system  evaluation  unnecessarily  hard.  A 
statistical  approach  reduces  the  information  to  a more 
understandable  form.  For  a given  magnification,  the  25 
spots  in  the  array  are  subjected  to  four  equal  increments, 
yielding  100  displacement  values  which  should  be  identical. 
Unfortunately,  the  displacements  returned  by  DPM  are  not  all 
equal.  For  the  sake  of  comparison,  the  standard  deviation 
of  the  100  displacement  values  is  calculated  and  normalized 
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by  the  accompanying  mean.  This  number  is  a measure  of  the 
amount  of  deviation  from  the  applied  constant  displacement 
and  serves  as  a vehicle  to  evaluate  spot  size. 

The  results  in  Table  1 show  that  there  is  a definite  lower 
limit  on  spot  size.  Although  difficult  to  say  precisely,  a 
spot  diameter  smaller  that  ten  pixels  is  considered  suspect. 


Table  1.  DPM  error  as  a function  of  spot  size. 


DIAMETER 

(PIXELS) 

STANDARD  DEVIATION 
DIVIDED  BY  THE  MEAN 

6 

0.019872 

9 

0.003052 

14 

0.003203 

16 

0.004553 

18 

0.004842 

24 

0.007174 

Ideally,  all  spots  in  the  array  should  give  the  same 
displacement  (for  rigid  body  motion) . System  and 
experimental  error  degrade  this  idealization  and  produce  a 
certain  amount  of  variance  from  spot  to  spot.  The  effect  of 
system  error  can  be  seen  by  examining  one  increment  of  the 
16  pixel  diameter  spots.  The  average  displacement  (over  the 
25  spots)  is  5.183  pixels  with  a standard  deviation  of 
0.0039  pixels.  Averaging  images  from  the  same  scene  reduces 
the  influence  of  random  system  noise  (cell  crosstalk, 
thermal,  electromagnetic,  etc.).  To  show  this,  four  images 
of  each  16  pixel  diameter  spot  scene  are  averaged  before 
using  DPM.  The  results  show  a significantly  reduced 
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standard  deviation  (0.00255  pixels)  which  suggests  a more 
accurate  measure  of  displacement  (5.181  pixels). 

Henceforth,  all  hybrid  experiments  use  four  image  averages. 

"False  system  displacement"  is  a result  of  the  variance 
in  displacement  between  spots  in  the  same  image  and  can  be 
evaluated  through  self-registration.  If  an  image  is 
registered  with  itself,  each  spot  should  measure  zero  rigid 
body  motion.  But  instead  the  system  displacements  are 
returned.  Table  2 reveals  the  displacement  errors 
introduced  by  the  system  and  represents  the  lower  accuracy 
limit  afforded  by  DPM  and  the  hardware. 


Table  2.  DPM  system  displacement  as  a function  of 
spot  size. 


DIAMETER 

(PIXELS) 

MEAN  SYSTEM 
DISPLACEMENT 
(PIXELS) 

6 

0.002251 

9 

0.001435 

14 

0.001063 

16 

0.001632 

18 

0.002117 

24 

0.001562 

The  above  errors  are  inherent  to  feature  extraction  and 
are  somewhat  abstract;  that  is,  no  definite  values  can  be 
used  to  discriminate  good  from  bad.  Search  region  over- 
expansion is  different  in  this  respect;  over-expansion 
defines  the  limit  to  the  displacement  increment  size.  For 
horizontal  rigid  body  motion,  one  expects  the  matching  to 


95 


fail  when  the  displacement  reaches  one  half  the  spot 
spacing.  This  limit  is  easily  verifiable  through  a simple 
experiment.  A 5x5  array  of  0.05  inch  diameter  spots  is 
again  tracked  through  0.025  (±0.001)  inch  increments;  but 
this  time  the  increments  are  continued  until  the  matching 
fails.  The  results  are  shown  in  Table  3 for  a spot  spacing 
of  0.2  inches  and  a spot  diameter  in  pixels  of  14.  As 
predicted,  the  matching  is  precise  until  a displacement 
greater  than  0.1  inch  is  applied.  The  limit  on  displacement 
increment  will  not  be  as  easily  identifiable  when  structural 
distortions  are  encountered,  but  a safe  rule  of  thumb  is  not 
to  exceed  one  fifth  of  the  smallest  spot  spacing. 

Table  3.  DPM  search  region  over  expansion. 


APPLIED 

MEASURED 

DISPLACEMENT 

DISPLACEMENT 

( inches) 

( inches) 

0.0 

0.00005 

0.025 

0.0249 

0.05 

0.0499 

0.075 

0.0751 

0.1 

0.0996 

0.125 

0.0286 

Displacement  pattern  matching  is  an  accurate,  efficient 
means  of  tracking  displacements.  All  DPM  calculations  are 
performed  on  an  80287-8  based  personal  computer  with  run 
times  averaging  about  1.75  minutes  per  increment.  For 
comparison,  DPM  installed  on  a VAX  11/750  completes  one 


96 


increment  in  less  than  12  seconds.  With  a few  exceptions, 
all  calculations  are  integer  greatly  enhancing  speed  and 
accuracy.  The  few  exceptions  are  single  precision  floating 
point  operations  in  the  centroid  and  standard  deviation 
calculations.  DPM  is  written  in  FORTRAN.  Although 
structured,  FORTRAN  has  notoriously  slow  input-output  (I/O) 
capabilities.  Ninety  percent  of  the  processing  time  is 
expended  reading  and  writing  image  information.  Experiments 
in  C (also  a structured  language)  show  significant 
reductions  in  I/O  time.  Reading  in  and  writing  out  an  image 
are  accomplished  in  7 seconds;  hence,  it  is  anticipated  that 
a C version  of  DPM  will  yield  personal  computer  run  times 
competitive  with  the  VAX. 

DPM  is  not  in  its  final  nor  optimum  form.  Improvements 
in  speed  can  be  achieved  through  the  use  of  different 
programming  languages,  such  as  C and  Assembler,  and 
different  hardware  [116].  Immediate  improvements  in 
accuracy  can  be  obtained  through  mathematical  camera 
corrections  [117]  and  system  calibration  [118].  However,  in 
its  present  stage  of  development,  displacement  pattern 
matching  is  certainly  comparable  to  other  sub-pixel 
registration  algorithms  in  accuracy;  and  it  frequently 
provides  a significant  speed  advantage  [90]. 


CHAPTER  IV 


HYBRID  DPM-EPBEM  TECHNIQUE 


Both  the  elastic-plastic  boundary  element  and 
displacement  pattern  matching  algorithms  have  been 
individually  verified,  but  the  ultimate  goal  here  is  to 
construct  a single  successful  stress  analysis  tool  by 
meshing  these  two  techniques.  Accordingly,  this  chapter 
presents  a hybrid  DPM-EPBEM  and  its  experimental 
verification.  The  following  sections  describe  the 
experimentation,  specimen  preparation,  and  verification 
tests. 


IV. 1 Working  Environment 

In  general,  combining  two  programs  written  in  the  same 
language  is  a trivial  matter.  However,  when  the  two 
programs  use  the  same  language  but  are  written  on  two 


97 


98 


):\JIffNflAP>del  scratch.dat 
:\JIff\HAP>del  cellook.dat 

:SJIffSMAP>pause  To  Exit  Graphics  Node  Type  * RESET'  . 
trike  a key  when  ready  . . . 
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PERCENT^  10.307 
THRESH0LD=118 


Figure  45.  Graphical  working  environment. 

different  computers  the  process  is  significantly  more 
complicated.  Here,  DPM  is  written  on  a PC's  Limited  286-8 
personal  computer  and  EPBEM  is  written  on  a Digital  VAX 
11/750.  These  two  computers  have  no  direct  means  of 
communicating  with  each  other,  so  a third  program,  PROCOMM, 
is  used  as  an  interpreter.  To  run  DPM-EPBEM,  the  image  of 
each  increment  is  digitized  and  the  resulting  digital 
representation  is  stored  in  such  a way  that  DPM  can  either 
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access  these  files  interactively  or  in  batch  mode.  In  any 
case,  the  user  controls  the  program  through  a graphical 
working  environment. 

The  working  environment  shown  in  Figure  45  provides 
critical  information  for  real  time  evaluation  of  the 
program's  progress.  A binary  pseudo-image  shows  the  scene 
currently  being  analyzed.  Included  in  this  picture  is  a 
window  encompassing  the  segmentable  portion  of  the  image. 

The  pseudo-image  is  produced  by  averaging  each  set  of 
sixteen  nearest  neighbors  and  then  thresholding.  The  result 
is  a 94H  x 121V  matrix  representation  of  a binary  image. 

This  representation,  as  opposed  to  a 384H  x 484V,  is 
necessitated  by  the  limited  dynamic  memory  afforded  by 
Microsoft's  DOS  and  the  WATFOR-87  FORTRAN  compiler.  On  the 
right  of  the  pseudo-image  is  a histogram  of  the  true  image 
with  the  location  of  the  threshold  indicated  by  the  vertical 
line.  The  corresponding  value  of  the  threshold  is  given 
below  the  histogram.  Also  provided  below  the  histogram  is 
the  percentage  of  pixels  with  the  most  frequent  gray  level. 
All  this  information  aids  in  evaluating  the  performance  of 
the  DPM. 

Above  the  pseudo  image  and  histogram  shown  in 
Figure  45,  is  a region  of  the  screen  reserved  for  any 
necessary  software  commands.  That  portion  of  the  screen  is 
called  the  editor  zone.  After  each  set  of  displacements  is 
determined  by  DPM,  the  communications  program,  PROCOMM,  is 
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evoked  through  a batch  file.  The  batch  file  logs  onto  the 
VAX,  transfers  the  displacement  files,  and  starts  EPBEM. 
During  this  process,  the  editor  zone  displays  each  batch 
command  as  it  is  executed.  Once  the  VAX  attains  control, 
all  of  the  graphics  are  temporarily  suppressed;  and  after 
EPBEM  has  completed  the  necessary  calculations,  the  VAX 
maintains  control.  VAX  control  enables  the  user  to  transfer 
the  stress  and  strain  information  back  to  the  host  computer 
(personal  computer) , to  return  to  the  host  computer,  or  to 
log  off.  Regardless  of  which  is  chosen,  all  of  the  desired 
stress  information  is  stored  on  the  VAX  and  can  easily  be 
retrieved  for  further  analysis. 


IV. 2 Experimental  Setup  and  Specimen  Preparation 

This  section  provides  a description  of  the  equipment 
and  of  the  general  experimental  setup.  Since  EPBEM  is  a 
two-dimensional  code,  only  plane  stress  or  plane  strain 
loading  conditions  can  be  examined.  All  specimens  are  of 
the  plane  stress  type  and  are  produced  from  a ductile 
aluminum  alloy  (1100-H14).  Loads  are  applied  by  a Material 
Testing  System  (MTS)  machine.  A schematic  of  the  setup  is 
given  in  Figure  46  and  a photograph  of  the  actual  system  in 
Figure  47.  The  testing  machine  is  manually  controlled  to 
produce  either  constant  load  or  constant  stroke. 
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Figure  46. 


Schematic  of  the  experimental  set  up. 
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As  described  in  Chapter  III,  DPM  depends  on  the 
application  of  high  contrast  spots  to  the  specimen.  The 
aluminum  specimens  are  painted  flat  white  in  order  to 
produce  the  high  contrast  background  for  the  black  spots. 


Figure  47.  Photograph  of  the  experimental  setup. 
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Initial  painting  attempts  resulted  in  the  paint  chipping  and 
cracking  under  load.  Since  DPM  perceives  chips  or  cracks 
as  spots,  one  must  undertake  means  to  stop  this  phenomenon 
from  occurring.  Each  specimen  is  first  thoroughly  sand 
blasted  to  remove  gross  deposits  of  aluminum  oxide.  They 
are  then  cleansed  in  a bath  of  isopropyl  alcohol.  Following 
air  drying,  each  specimen  is  scrubbed  carefully  with 
AMCHEM's  Alumiprep  33  (stock  no.  DX533)  to  remove  the 
remaining  corrosive  oxidation  and  to  chemically  etch  the 
surface.  The  final  treatment  is  a generous  scrub  and  bath 
in  AMCHEM's  Alodine  1201  (stock  no.  DX503).  This  treatment 
chemically  stabilizes  the  aluminum  surface  and  provides  a 
long  lasting  paint  adhesive.  An  additional  advantage  of 
this  procedure  is  that  it  leaves  a distinctive  amber  stain 
where  the  coating  process  is  successful.  The  specimens  are 
then  air  dried  and  painted  with  flat  white  KRYLON  enamel. 

The  entire  process,  from  start  to  finish,  reguires  about  ten 
minutes  per  specimen  and  provides  a painted  surface  which 
does  not  chip  or  crack  until  near  catastrophic  failure  loads 
are  applied. 

Spot  application  followed  an  iterative  process  until 
finding  a satisfactory  procedure.  The  primary  concerns  here 
are  uniformity  in  spot  diameter  and  darkness.  Best  results 
are  obtained  by  using  typical  "rub-ons"  commonly  found  in 
office  supply  stores;  the  grammatical  periods  are  far  more 
consistent  in  area  and  darkness  than 
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any  attempted  painting  or  staining  procedures.  The  periods 
are  simply  rubbed  on  at  the  locations  where  the  displacement 
is  desired. 

Three  plane  stress  experiments  are  presented  in  order 
of  increasingly  complex  states  of  stress.  The  first  test 
examined  is  a thin  sheet  under  uniaxial  tension,  the  second 
test  is  a thin  rectangular  sheet  with  a central  hole, 
subjected  to  tensile  end  loads  (perforated  strip  tensile 
test)  and  the  final  test  is  a thin  rectangular  sheet  with 
symmetric  ninety  degree  V-notches.  The  1100-H14  aluminum 
used  for  all  specimens  has  the  following  material 
properties:  E=10,713.  kpsi,  j/  = .33,  and  Yo=15,000  psi. 

EPBEM  reguires  the  uniaxial  stress  strain  curve  in  the 
iterative  solution.  To  accommodate  EPBEM,  the  uniaxial  true 
stress-engineering  strain  curve  is  found  experimentally  by 
applying  known  loads  to  a uniaxial  test  specimen  with  the 
MTS  machine  and  monitoring  a one  inch  gage  length  with  the 
MTS  supplied  extensometer . The  stress-strain  data  is  fitted 
with  two  piecewise  continuous  curves.  The  elastic  and  knee 
portion  are  fitted  with  the  Ramberg-Osgood  Law  [118];  and 
the  work  hardening  portion  is  approximated  with  a linear 
fit.  The  approximate  and  experimental  stress-strain  curves 
are  given  in  Figure  48. 
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Figure  48.  Uniaxial  stress-strain  curve  for  1100-H14 
aluminum. 


IV. 3 Calibration 


Prior  to  each  experiment,  the  digitizing  system  is 
calibrated  carefully  so  that  DPM  results  are  accurately 
converted  from  pixels  to  units  of  length.  The  calibration 
constant  is  a function  of  magnification  so  that  the  first 
step  in  calibrating  the  system  is  to  set  the  desired 
magnification  and  to  keep  it  fixed  for  all  load  increments. 
Next,  a 5x5  array  of  black  dots  is  placed  in  the  focal  plane 
and  given  five  equal  vertical  displacements;  the  array  is 
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digitized  after  each  step.  Since  all  of  the  displacements 
are  equal,  a linear  relation  exists  between  the 
displacements  returned  by  DPM  (in  pixels)  and  the  measured 
displacements;  the  slope  of  this  curve  is  the  calibration 
constant.  The  steps  to  calculating  the  calibration  constant 
are  as  follows: 

1)  Apply  each  load  increment. 

2)  Store  the  four  image  averages  of  each 
scene. 

3)  Calculate  the  average  rigid  body 
displacement  of  the  25  spots  in  each 
scene . 

4)  Plot  the  average  displacement  (in  pixels) 
versus  the  measured  displacement  (in  units 
of  length) . 

5)  Use  a linear  least  squares  fit  to  find  the 
slope. 

The  imaging  array  of  the  digitizing  board  used  in  these 
experiments  is  rectangular;  therefore,  the  y-calibration 
constant  differs  from  the  x-calibration  constant.  To 
accommodate  the  rectangular  imaging  array,  the  above 
calibration  procedure  is  repeated  with  horizontal 
displacements.  A typical  experiment  produced  x-  and  y- 
calibration  constants  of  369.5  and  523.3  pixels/inch, 
respectively.  With  system  resolution  on  the  order  of 
.1  pixels,  the  displacement  resolution  is  on  the  order  of 
.0001  inch. 
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IV. 4 Uniaxial  Tension  Test 


Before  describing  the  hybrid  experiments,  it  should  be 
noted  that  no  attempt  is  made  to  relax  machining  induced 
work  hardening  in  the  specimens.  Although  care  is  taken  to 
use  sharp  tools  and  slow  cutting  speeds,  no  annealing  is 
performed.  The  specimens  for  all  experiments  are  produced 
from  the  same  sample  of  aluminum;  and  the  uniaxial  stress 
strain  curve  for  each  specimen  is  assumed  to  be  the  same. 

From  the  least  complex  test  of  DPM-EPBEM,  uniaxial 
tension,  one  expects  to  uncover  any  limitations  of  the 
hybrid  technique.  With  this  in  mind,  the  hybrid  analysis  of 
a uniaxial  tension  specimen  is  described.  The  dimensions  of 
this  plane  stress  specimen  are  nine  inches  tall,  one  inch 
wide  and  .0625  inches  thick  (all  specimens  examined  have 
these  nominal  dimensions) . The  spots  are  applied  to  the 
specimen  as  indicated  in  Figure  49  and  as  shown  in 
Figure  50.  Two  interior  triangular  cells  are  used  to  model 
the  region  enclosed  by  the  spots,  with  DPM  providing 
displacement  boundary  conditions  at  all  but  the  free 
surface.  One  should  note  that  only  the  region  encompassed 
by  DPM  is  modelled  by  the  boundary  element  method  and  that 
the  size  and  location  of  this  region  are  completely 
arbitrary.  Figure  39  shows  a plot  of  the  x-component  of 
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Figure  49.  Schematic  of  the  spot  locations  on  the 
uniaxial  tension  specimen. 

displacement  which  is  calculated  by  DPM  at  (x=0.25,y=0. ) . 

The  results  conform  to  the  shape  one  expects  for  a load- 
displacement  curve  of  weakly  work  hardening  materials.  A 
little  more  convincing  is  the  plot  shown  in  Figure  52,  which 
provides  the  stress-strain  curve  calculated  by  DPM-EPBEM 
along  with  the  experimental  version.  Two  interesting 
characteristics  of  note  are  the  good  correlation  between  the 
MTS  and  DPM-EPBEM  results  and  the  limited  useful  range  of 


DPM-EPBEM. 
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Figure  50.  Photograph  of  the  spots  applied  to  the 
uniaxial  test  specimen. 

On  the  lower  end,  the  DPM-EPBEM  curve  is  unreliable 
until  approximately  90%  of  the  proportional  limit.  This  is 
a result  of  the  rigid  body  motion  introduced  by  imperfect 
alignment  of  the  specimen  in  the  MTS  wedge  grips.  Although 
the  specimens  are  aligned  as  carefully  as  possible,  no 
aligning  gages  are  used.  As  the  initial  load  steps  are 
applied,  the  specimen  adjusts  itself  to  evenly  distribute 
the  load.  The  redistribution  introduces  both  rigid  body 
rotation  and  translation.  A second  source  of  rigid  body 
motion  comes  from  the  specimen  slipping  while  the  wedge 
grips  are  locking  down.  This  motion  also  occurs  during  the 
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Figure  51.  Load-horizontal  displacement  curve  for 
(x=0 , y=. 25) . 
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Figure  52 


Experimental  and  DPM-EPBEM  uniaxial  stress 
strain  curve. 
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initial  loading  steps  and  although  it  is  less  pronounced 
than  the  rigid  body  rotation  from  specimen  misalignment, 
this  slippage  is  certainly  significant. 

DPM  is  able  to  discern  the  rigid  body  motion  from  both 
sources  and  passes  this  rigid  body  motion  information  to 
EPBEM.  Unfortunately,  EPBEM  expects  the  constrained 
displacements  and  has  no  means  to  separate  out  the  rigid 
body  motion.  Therefore,  the  boundary  element  routine 
calculates  erroneous  stress  and  strain  values. 

Whereas  DPM  is  responsible  for  the  lower  performance 
limit,  EPBEM  is  responsible  for  the  upper  performance  limit. 
EPBEM  expects  a work  hardening  material  and  has  difficulty 
converging  to  a solution  for  a perfectly  plastic  material. 
The  upper  limit  results  from  the  fact  that  1100-H14  aluminum 
behaves  in  a almost  perfectly  plastic  manner.  At  high  a 
stress,  the  flat  portion  of  the  stress  strain  curve  causes 
the  iterative  solution  procedure  to  become  unstable; 
essentially,  the  initial  strain  procedure  fails  to  converge 
to  a single  value.  It  is  believed  that  the  load  increments 
in  the  flat  portion  of  the  stress-strain  curve  produce 
strain  increments  too  large  for  the  initial  strain  solution 
procedure  to  assimilate.  An  upper  performance  limit  is 
typical  in  all  of  the  experiments  performed  but  it  varies 
with  loading  conditions  and  geometry.  Still,  the  DPM-EPBEM 
upper  performance  limit  usually  occurs  at  an  equivalent 
stress  of  16,500  psi.  Overall,  the  uniaxial  tension  test 
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exhibits  dependable  results  within  a specific  equivalent 
stress  range  and,  as  anticipated,  unmasks  basic  weaknesses 
of  the  DPM-EPBEM  hybrid  technique. 


IV. 5 Perforated  Strip  Tensile  Test 

The  uniaxial  tensile  test  produces  comforting  results. 
In  addition,  the  results  establish  both  the  accuracy  and  the 
performance  range  of  the  hybrid  technique.  This  next  test 
seeks  the  solution  to  a more  realistic  engineering  problem 
and,  among  other  important  details,  the  perforated  strip 
exposes  DPM-EPBEM  to  a stress  state  containing  high  stress 
gradients.  The  stresses  calculated  by  the  hybrid  technique 
are  compared  to  both  ANSYS  finite  element  and  Theocaris  and 
Marketos'  experimental  results  [120]. 

The  spots  are  applied  to  the  specimen  at  the  locations 
specified  in  Figure  53  and  shown  in  Figure  54;  where  the 
domain  enclosed  by  the  spots  is  modelled  as  seen  in  Figure 
55.  A careful  comparison  between  Figure  55  and  Figure  53 
reveals  that  not  all  of  the  DPM  spots  and  EPBEM  nodes  are 
found  at  the  same  locations.  The  finite  size  of  the  spots 
establishes  a practical  limit  to  the  spot  spacing; 
therefore,  DPM-EPBEM  uses  quadratic  interpolation  to  find 
the  displacements  at  nodes  which  do  not  coincide  with  spots. 
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Figure  53.  Schematic  of  the  spot  locations  on  the 
perforated  strip. 
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Figure  54.  Photograph  of  the  spots  applied  to  the 
perforated  strip. 

The  homogeneous,  isotropic  hardening  option  (von  Mises 
yield  condition)  and  two-dimensional,  four  point, 
isoparametric  plate  elements  (two  degrees  of  freedom)  are 
used  in  the  finite  element  solution  to  this  problem.  ANSYS 
follows  an  initial  stress  [121]  solution  procedure, 
modelling  the  stress-strain  curve  with  five  piecewise  linear 
segments.  One  symmetric  quarter  of  the  specimen  is  modelled 
(Figure  56)  and  is  loaded  with  50  psi  traction  increments, 
added  to  the  initial  traction  of  3000  psi.  The  load  steps 
are  applied  until  a value  of  6300  psi  is  reached.  The  DPM- 
EPBEM,  ANSYS,  and  Theocaris ' stress  concentration  factors  at 
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Figure  55. 


Boundary  element  domain  discretization  of 
the  perforated  strip. 
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Figure  56.  Finite  element  grid  for  the  perforated 
strip. 

(x=0.25,  y=0.0)  as  a function  of  load  step  are  given  in 
Figure  57.  From  here  on,  the  stress  concentration  factor 
(SCF)  is  defined  as  the  current  y-component  of  stress 
divided  by  the  applied  stress.  The  finite  element  and  DPM- 
EPBEM  results  compare  well,  with  DPM-EPBEM  exhibiting  one 
larger  excursion  at  the  sixth  load  step.  No  immediate 
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Figure  57.  Plot  of  stress  concentration  factor 
versus  load  step  for  the  perforated 
strip. 


explanation  for  this  anomaly  is  apparent.  Since  the  next 
increment's  result  does  not  seem  to  be  affected,  one  might 
presume  that  some  grip  slippage  had  occurred. 

Perhaps  even  more  comforting  are  the  experimental 
points  produced  by  Theocaris  and  Marketos  [120].  The 
material  they  tested  was  also  an  aluminum  alloy.  It 
exhibited  a similar  uniaxial  stress-strain  curve  to  Figure 
51  with  a higher  yield  strength  (34,500  psi) . Although  the 
SCF  curves  are  not  identical,  which  is  not  to  be  expected, 
the  purely  experimental,  hybrid  and  purely  numerical 
behaviors  are  remarkably  similar. 
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The  difference  in  the  response  could  be  caused  by  one 
of  many  reasons.  Theocaris  and  Marketos  use  a bi-linear 
approximation  of  the  uniaxial  stress-strain  curve  in  their 
analysis  whereas  the  DPM-EPBEM  uses  a combination  of 
Ramberg-Osgood  and  linear  fits.  Theocaris  and  Marketos  make 
no  mention  of  whether  they  did  or  did  not  anneal  their 
specimens,  nor  do  they  discuss  their  loading  apparatus  in 
any  detail.  Therefore,  one  would  not  necessarily  expect 
perfect  agreement. 

Other  interesting  information  produced  by  the  hybrid 
technigue  is  the  load-displacement  curve  for  (x=0,y=.25) 
which  is  given  in  Figure  58.  The  load-displacement  results 
mimic  (in  shape)  the  load-maximum  strain  curves  reported  by 
Zienkiewicz  [122].  Although  not  definitive,  they  certainly 
are  supportive  of  DPM-EPBEM  reliability.  In  his  study, 
Zienkiewicz  used  the  material  properties  reported  by 
Theocaris  and  Marketos. 
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Figure  58. 


Load-displacement  plot  for  (x=0.,y=.25)  on 
the  perforated  strip. 
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Figure  59 
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A final  plot.  Figure  59,  demonstrates  the  extent  of  the 
plastic  zone  by  giving  the  SCF  at  various  radial  locations 
from  the  root  of  the  hole  (for  ty=6000  psi) . The  slight 
inflections  at  x=0.4  and  x=0.275  are  apparent  in  all  of  the 
graphs  provided  by  Theocaris  and  Marketos  [120]  and  are 
therefore  reassuring. 


IV. 6 V-notched  Specimen  Tensile  Test 

The  correlation  between  DPM-EPBEM,  finite  element 
analysis  and  the  purely  experimental  results  establish  a 
footing  from  which  a more  complex  problem  can  be  approached. 
The  concluding  experiment  attempts  to  subject  the  hybrid 
technique  to  a stress  state  truly  representative  of  common 
engineering  problems,  a notched  tensile  specimen  [123]. 

This  geometry  (Figure  60)  produces  high  stress  gradients 
localized  near  the  root  of  the  notch  with  a complex  stress 
field  in  the  surrounding  material  [124].  It  should  be  noted 
that  there  is  a small  machine-induced  rounding  of  the  notch 
root;  therefore,  the  notch  can  not  be  considered  sharp  in 
the  ideal  sense.  Once  again,  the  ANSYS  finite  element  code 
and  selected  published  results  are  compared  to  those 
produced  by  the  DPM-EPBEM  hybrid  technique. 
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Figure  60.  Schematic  of  the  spot  locations  on  the 
V-notched  tensile  specimen. 
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As  before,  the  spots  are  positioned  according  to 
Figure  60  and  displayed  in  Figure  61.  A sketch  of  the 
boundary  element  grid  can  be  found  in  Figure  62  and  one  of 
the  finite  element  grid  in  Figure  63.  Again,  the 
homogeneous,  isotropic  hardening  option  and  two- 
dimensional,  four  point  isoparametric  elements  are  used  in 
the  finite  element  solution  of  this  problem. 


Figure  61.  Photograph  of  the  spots  applied  to  the 
V-notched  tensile  specimen. 
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Figure  62.  Boundary  element  domain  dicretization  of 
the  V-notched  tensile  specimen. 

Both  the  SCF  vs  load  increment  (Figure  64)  and  the  load 
-displacement  curve  for  (x=0,y=.47)  (Figure  66)  are  provided 
for  evaluation  of  the  performance  of  DPM-EPBEM . The  finite 
element  stress  concentration  factor  plot  is  calculated  by 
first  applying  an  end  load  of  2000  psi,  and  by  then  adding 
load  increments  of  50  psi  until  reaching  a value 
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Figure  63.  Finite  element  mesh  for  the  grid  for  the 
V-notched  tensile  specimen. 

of  5400  psi.  In  Figure  64,  the  hybrid  and  FEM  results 
differ  markedly.  The  first  point  of  deviation  is  at  first 
yield.  DPM-EPBEM  predicts  that  the  first  yield  will  occur 
at  an  applied  traction  of  2800  psi,  whereas  ANSYS  predicts 
the  first  yield  at  3100  psi.  One  can  not  say  which  is  a 
better  solution  but  domain  discretization  does  have  an 
impact  on  the  comparison.  If  the  notch  is  perfectly  sharp, 
then  plastic  deformation  initiates  when  the  first  load  is 
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Figure  64.  Plot  of  stress  concentration  factor  versus 
load  step  for  the  V-notched  tensile 
specimen. 

applied.  The  deformation  field  is  so  localized  that  many 
internal  cells  placed  near  the  notch  tip  would  be  needed  to 
discern  the  motion.  The  stress  calculated  at  each  cell  is 
an  average  representation.  Therefore,  at  the  notch  root, 
the  larger  cells  produce  smaller  average  values.  A 
comparison  shows  that  the  finite  element  grid  has  three  less 
cells  near  the  notch  root  than  does  the  boundary  element 
grid;  therefore,  it  is  reasonable  to  expect  that  the  first 
yield  predictions  by  BEM  and  FEM  will  be  different. 
Experimental  error  also  contributes  the  different  initial 
yield  predictions  of  DPM-EPBEM  and  FEM.  The  erratic  hybrid 
SCF  curve  suggests  that  some  of  the  DPM  displacements  may  be 
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Figure  65.  Photograph  of  the  V-notched  tensile 

specimen  before  loading  and  after  failure. 
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in  error.  If  this  is  so,  then  the  early  increments  would 
equally  be  in  error.  Regardless  of  the  reasons  for  the 
difference  in  the  predictions,  the  qualitative  similarity  is 
apparent  with  the  quantitative  results  differing  at  most  by 
14.9%. 

Other  important  considerations  deal  primarily  with  the 
difference  between  the  ideal  loading  conditions  in  the 
numerical  simulation  and  the  actual  loading  conditions  in 
the  hybrid  technique.  Examination  of  the  post  yield 
specimen  (Figure  65)  indicates  that  yield  did  not  occur 
simultaneously  at  both  notch  roots.  This,  in  turn,  implies 
that  the  specimen  was  not  perfectly  aligned  in  the  MTS 
grips,  precipitating  a non-symmetric  stress  state.  In  such 
a case  as  this,  one  cannot  expect  the  FEM  solution  to 
duplicate  the  experimental  results. 

Although  no  comprehensive  experimental  investigation  of 
the  V-notched  test  specimen  is  available,  inferences  from 
results  presented  by  several  previous  experimental  and 
numerical  studies  lend  credence  to  the  results  presented 
here.  The  load  displacement  curve  given  in  Figure  66 
qualitatively  compares  well  with  a similar  numerical  curve 
reported  by  Telles  [18]  (page  171)  and  with  the  experimental 
curve  reported  by  Findley  and  Drucker  [125]  (1/16  inch 

specimens) . In  both  cases,  the  material  is  an  aluminum 
alloy  which  is  different  from  the  material  used  here. 
Additionally,  the  aspect  ratios  and  the  net  cross  sectional 
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Figure  66.  Load-displacement  plot  for  (x=0,y=.46)  on 
the  V-notched  tensile  specimen. 
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Figure  67.  Plot  of  stress  concentration  factor  versus 
x-position  on  the  V-notched  tensile 
specimen. 
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areas  of  their  specimens  differ  from  those  used  in  the 
present  experiments.  Taking  these  facts  into  consideration 
and  the  fact  that  EPBEM  reached  its  upper  performance  limit 
before  the  plastic  limit  load.  Figure  66  exhibits  the 
characteristic  linear  load-displacement  curve  in  the  pre- 
limit load  regime.  Linear  behavior  from  a plastically 
deforming  body  may  seem  contradictory,  but  in  this  case  it 
is  reasonable.  During  the  early  stages  of  plastic 
deformation  (as  is  the  case  with  DPM-EPBEM)  the  plastic  zone 
is  localized  in  the  region  near  the  root  of  the  notch.  The 
surrounding  elastic  material  confines  the  plastic  zone  and 
effectively  stems  any  large  plastic  deformation.  The 
regions  away  from  the  root  are  essentially  unaware  of  the 
plastic  flow  until  the  loads  are  high  enough  to  produce 
substantial  redistribution  of  the  internal  stress  state  of 
the  specimen  and,  in  turn,  a larger  plastic  flow  [125]. 

Finally,  the  stress  concentration  factor  along  the  y=0 
axis  (Figure  67)  follows  the  expected  course  by  increasing 
as  one  moves  closer  to  the  root.  The  EPBEM  upper 
performance  limit  (5400  psi)  prevents  a fully  developed 
plastic  zone  and  hence  a flattening  out  of  the  SCF  curve  is 
not  present.  The  SCF  versus  position  plot  is  taken  for  an 
applied  load  of  4700  psi,  a load  at  which  the  plastic  zone 
is  still  confined  to  the  notch  root. 


CHAPTER  V 


CONCLUSIONS 


The  two-dimensional  numerical-experimental  hybrid 
technique  described  herein  combines  context-free  syntactic 
pattern  recognition  (DPM)  and  elastic-plastic  boundary 
element  methods  (EPBEM)  to  produce  a useful,  non-destructive 
stress  analysis  tool.  Within  a certain  load  range,  all  of 
the  presented  DPM-EPBEM  results  compare  well  with  finite 
element  and  purely  experimental  (where  available)  solutions. 
The  DPM-EPBEM  solution  for  the  V-notched  tensile  specimen 
exhibits  a more  erratic  response  than  the  responses  of  the 
other  two  experiments.  Since  the  V-notched  tensile  specimen 
represents  a state  of  stress  that  is  as  complex  as  one  would 
expect  to  encounter  in  engineering  problems,  an 
understanding  of  the  results  of  that  test  is  paramount  to 
evaluating  DPM-EPBEM' s usefulness  as  a stress  analysis  tool. 

DPM-EPBEM  works  well  in  load  ranges  where  rigid  body 
motion  is  negligible  when  compared  to  the  displacements 
produced  by  distortion.  Rigid  body  motion  in  the  initial 
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stages  of  loading  (approximately  90%  of  the  proportional 
limit)  degrades  the  accuracy  of  the  DPM-EPBEM  hybrid 
technique.  This  lower  accuracy  limit  can  be  improved  by 
incorporating  rigid  body  motion  compensation  in  DPM  or 
EPBEM. 

In  instances  when  DPM-EPBEM  returns  questionable 
results,  a possible  question  might  be  raised  about  how  to 
guarantee  that  the  applied  spots  are  in  the  exact  locations 
as  specified  in  Figures  49,  53  and  60.  The  accuracy  of  the 
spot  application  technique  described  in  Chapter  IV  is 
limited  by  human  error.  In  most  cases,  the  accuracy  of  the 
spot  location  technique  is  approximately  ±.001  inch. 
Surprising  as  it  may  seem,  mirroring  the  schematics  is  not 
that  critical.  The  boundary  element  code  does  receive 
displacement  information  from  DPM,  but  it  also  receives  and 
stores  the  locations  of  the  spot  centroids  in  the  undeformed 
image.  Using  this  information,  DPM-EPBEM  interpolates  (or 
extrapolates  as  the  case  may  be)  the  displacement 
information  to  the  BEM  node  locations.  This  is  not  to  say 
that  errors  are  not  introduced  by  misalignment  of  the  spots. 
If  the  spots  are  to  be  located  along  a horizontal  line,  then 
misalignment  in  the  vertical  direction  can  play  a 
significant  role.  However,  DPM-EPBEM  is  very  forgiving  with 
respect  to  spot  placement  errors. 
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Given  that  the  spots  are  perfectly  aligned,  the  ability 
of  DPM  to  measure  the  displacements  accurately  is  certainly 
paramount  to  the  success  of  the  hybrid  technique.  In  the 
images  used  in  this  investigation,  typical  values  for 
signal-to-noise  ratio,  gray  level  difference  and 
displacement  resolution  are  200,  100  and  0.0001  inch, 
respectively.  In  many  situations  there  are  regions  in  the 
tested  specimens  where  this  resolution  is  insufficient.  As 
an  example,  the  stress  concentration  factors  for  the 
perforated  strip  and  the  V-notched  tensile  specimens  are 
calculated  at  the  boundary  of  the  respective  cutout  roots. 
Therefore,  based  on  the  parametric  study  in  section  II. 5, 
the  DPM  displacement  error  plays  an  important  role  in  the 
SCF  results.  Both  cases  include  localized  plastic  zones  and 
the  same  nominal  dimensions.  Because  of  these  localized 
regions,  one  anticipates  that  the  displacement  in  the  region 
near  the  cutout  roots  governs  the  accuracy  of  the  DPM 
displacements.  Minor  displacement  errors  away  from  the 
cutout  roots  do  not  have  the  influence  of  displacement 
errors  in  the  plastic  zone.  The  plots  in  Figures  57  and  64 
follow  the  above  reasoning;  if  the  resolution  for  both 
experiments  is  the  same,  then  only  regions  of  small 
displacements  (on  the  order  of  the  resolution  limit)  can 
account  for  the  difference  in  the  behavior  of  the  SCF  plots. 
This  conclusion  immediately  draws  attention  to  the 
displacement  in  the  regions  surrounding  the  cutout  roots. 
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The  plastic  flow  in  the  perforated  strip  is  not  as  confined 
as  it  is  in  the  V-notch  specimen.  Therefore,  the 
displacement  at  the  root  of  the  circular  cutout  experiences 
larger  and  more  measurable  displacements.  The  displacements 
near  the  V-notch  are  much  more  confined  by  the  surrounding 
elastic  region  which  holds  the  displacements  nearer  to  the 
resolution  limit  of  DPM.  As  a result,  the  SCF  plot  for  the 
V-notch  specimen  is  much  more  erratic. 

One  should  note  that  the  specimens  used  in  the 
evaluation  of  DPM-EPBEM  are  smaller  (nominal  dimensions  are 
one  inch  wide,  nine  inches  tall  and  .0625  inch  thick)  than 
one  would  expect  to  encounter  in  many  true  engineering 
problems.  The  displacements  in  larger  specimens  will  be 
bigger  and  more  measurable  than  the  displacements 
encountered  here.  This  argument  holds  for  localized  plastic 
zones;  they  too  will  encompass  larger  areas,  providing  more 
sizable  displacements.  Therefore,  in  larger  structures,  the 
DPM-EPBEM  is  less  likely  to  encounter  the  resolution  limit 
of  DPM.  Avenues  for  improving  the  DPM  resolution  are 
available  for  problems  dictating  measurements  on  the  order 
of  the  resolution  limit.  The  interpolation  of  the 
displacement  information  from  spot  locations  to  node 
locations  has  the  effect  of  smoothing  out  the  influence  of 
the  poor  data  near  the  roots  of  the  cutouts.  This  suggest 
that  smoothing  technigues  tailored  to  these  types  of 
problems  may  be  useful.  The  real  solution,  however,  is  to 
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improve  the  resolution  of  DPM.  This  can  be  accomplished  by 
increasing  the  calibration  constants  through  magnification 
or  with  better  hardware,  or  by  compensation  for  system 
errors . 

Another  area  which  needs  improvement  deals  with  the 
elastic-plastic  boundary  element  method.  The  convergence  of 
the  initial  strain  elastic-plastic  boundary  element  solution 
is  sensitive  to  load  increment  step  size.  As  seen  in  the 
experiments  presented  in  Chapter  IV,  it  is  not  always 
possible  to  apply  the  loads  in  such  way  as  to  stay  within 
the  realm  of  incremental  plasticity  and  thereby  guarantee 
convergence.  Therefore  a solution  procedure  which  is  less 
sensitive  to  increment  size  would  be  more  desirable. 
Fortunately,  the  initial  stress  solution  approach  has  just 
this  quality  [126].  It  is  anticipated  that  an  initial 
stress  solution  procedure  would  greatly  enhance  the  load 
range  in  which  this  hybrid  technique  is  applicable. 

The  honest  and  somewhat  critical  evaluation  of 
DPM-EPBEM  is  given  in  the  above  paragraphs  so  that  the 
limitations  of  this  technique  can  be  understood  and  so  areas 
of  further  research  can  be  pinpointed.  But  even  in  its 
present  stage  of  development,  the  hybrid  experimental- 
numerical  procedure  described  in  this  report  is  a useful 
stress  analysis  tool.  The  good  results  presented  in 
Chapter  IV  are  one  indication  of  the  promise  that  the  DPM- 
EPBEM  procedure  holds.  A second  indication  of  the  promise 
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of  this  technique  is  that  it  is  fully  automated.  Finally, 
the  equipment  required  to  produce  a hybrid  method  like  the 
one  described  is  standard  at  most  engineering  and  research 
facilities.  By  incorporating  the  above  suggestions,  the 
displacement  pattern  matching,  elastic-plastic  boundary 
element  hybrid  method  can  be  shaped  to  provide  good 
solutions  over  an  even  wider  load  range. 


APPENDIX 


The  two  volume  integrals  encountered  in  elasto-plastic 
boundary  element  methods  appear  in  the  boundary  constraint 
equation  (Equation  11-15)  and  the  interior  strain  equation 
(Equation  11-40) . Here  these  integrals  are  identified  and 
their  semi-numerical  integration  is  described. 

The  constraint  equation  volume  integral  is  written  as 


where 

ASijk  ~ sijk(r)^v 

V 

sijk(r)  ~ ^ [c2(5ijr,k  + 5ikr,j  " 5jkr,i)  + r,ir,jr,k]  (A  1) 


with 

Cl  — 1 

4?r  ( l—i/ ) 

c2  = (l-i/)  , 

rj_  = - £ i and 

r/i 

axi 
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To  change  from  plane  strain  to  plane  stress  replace 
v by  . As  described  in  Chapter  II,  the  volume 

integral  is  evaluated  over  each  triangle  (Figure  8 and  9) 
and  algebraically  summed  to  find  the  integral  over  each 
interior  cell.  It  is  instructive  to  show  this  integration 


over  a representative  triangle  so  as  to  understand  its 


Figure  68.  Coordinate  system  for  the  integration  over 
an  arbitrary  triangular  cell. 
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A transformation  to  cylindrical  coordinates  simplifies 

I 

the  integration. 


AS 


ijk  (r)  = lim 


'0b 


0; 


fR(0) 


sijk  (r,<£)rdrd <f> 


(A-2) 


where  r,  <f> , 9a,  <j> b,  R (<f>)  and  Sijk  are  defined  in  Figure  68. 
The  integration  is  carried  out  with  respect  the  line 
segment,  D,  perpendicular  to  R(<£).  This  introduces  the  new 
coordinate  system  r 1/  2 • BY  defining  the  derivatives  of  r 

in  cylindrical  coordinates,  Sfjk  is  properly  recast. 


ar  a_r  d£_i  , dr  ar 
+ 2 
axi  ari  axi  ar2  3xi 


(A-3) 


As  seen  from  Figure  68,  — r = cos <f> , — = siruf> , = e^i 

ari  ar 2 axi 

and  = e2i;  where  eij  are  the  direction  cosines  between 

axi 

the  x and  r coordinate  systems.  Therefore, 

r,i  = cos <f>  e^i  + sin<£  e2i,  (A-4) 

and  combining  Equations  (A-l) , (A-2)  and  (A-4)  gives 


ASijk  - CiC2 


+ C!C2 


‘0b 

R(4>  ) a ij  (cos0  elk  + sin0e2k)d0 
0a 

•0b 

R(0  ) a ik(cos5  eij  + sin  <t>  e2j)d^ 

J 0a 


‘0b 

R(0)  5 jk(cos^eli  + sin0  e2i)d<£ 
0a 


(A-5) 


- CiC2 
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+ 2C]_ 


‘0b 

[R(<f>)  (cos^exi  + sin<^e2i)  (cos^exj 
0a 


+ sin^e2j) 


+ (cos0eik  + sin^  e2k)]d<£. 


From  the  geometry  shown  in  Figure  68,  <f>  = 0+a , d<f>=d9 
and  R{</>)  = D/cos^. , and  therefore  the  integral,  ^ijk  is 


^ijk  - dc1c2  ( (elk5  ij  + eljsik  “ elj5jk) 

+ tan0(e2kSij  + e2jfiik  - e2i§jk)} 

+ DC2{2cos25e1je1jelk 

+ 2cos^sin0 (e2ie1jelk  + e^e^e^  + e1je1je2k) 

+2sin20 (e2ielje2k  + elie2je2k  + e2ie2jelk)  (A_6) 

+ 2sin20 tan0 e2ie2 je2k} . 

. * 

Upon  integration  Equation  (A-5)  reduces  the  integral  of  Sijk 
over  an  arbitrary  triangle. 


A2ijk  = c1c2d(  (elk5  ij  + elj5ik  “ eli5jk)Il 

+ (e2k5ij  + e2jfiik  ~ e2isjk)I2}  + DC2elielj elkI3 
+ DC2 (e2ieijelk  + eije2jelk  + e1ie1je2k)I4 
+ DC2(e2ielje2k  + eiie2je2k  + e2ie2jelk)I5 
+ DC2e2ie2je2kI6,  (A~7) 
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where 


3 

II 

H 

H 

"0b 

d0  = #b  “ O a ' 
0 a 

12  = 

tan0d0  = -ln(  b)  _ in(  b)  , 

6a  cos0 a r2 

13  = 

‘0b  , 

2cos20d0  = (^b“0a)  + “(sin2^b-sin20 a) , 

0a  2 

14  = 

‘0  b 

2cos0sin0d0  = sin20b_sin20 a > 
0 a 

15  = 

■0b 

2sin20d0  = ( 6 b-0  a)  + ~(sin20b“sin20 a) 
0a  2 

16  = 

‘0b 

2tan0Sin20d0  = 21n(— b)  + cos20 b-cos20 a 
0 a ra 

The  volume  integral  associated  with  the  total  strain 
equation  is  simply  the  derivative  of  Equation  (A-7)  (see  the 
development  of  Equation  11-40) , that  is 

AE^ijk  = 7T  (A2ijk)  • (A-8) 

a 

Here  the  derivative  is  taken  with  respect  to  the  load  point. 
If  £ is  perturbed  then  D and  its  associated  angles  change; 
Ci,  C2  and  e^j  remain  fixed.  Evaluating  Equation  (A-8)  is 
now  reduced  to  finding 
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I[n  = (DIi)  where  i = 1,  2,  ...  6, 

30 


which  by  the  chain  rule  expands  to 


•ii  = 


3D 

aO 


li  + 3Ii 


^ D 


ae  a £ $ 


(A-9) 


The  first  term  — = -ei o and  — in  the  second  term  is 

30  3 O 


evaluated  as  follows: 

6 = tan-1  (^2)  - a , 

O 


M.  _ _ £2 

ae  R ' 


(A-10 ) 


but  i2  ~ Rsin(i9+oc)  , after  expanding  and  substituting  into 
Equation  (A-10) 


a 6 sing  cos a - cosflsina 


3£. 


R 


similarly 


a 6 cosgcosa  - smasing 


a?. 


R 


By  making  use  of  the  direction  cosines 


dj_ 

aO 


e2£COS0  ~ eii5in 9 


R 
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Now,  by  substituting  the  above  results  into  Equation  (A-9) 
AElijk  can  be  written  for  a representative  triangle 

AEiijk  = clc2(elk5ij  + elj5ik  “ eli5jk)Il 
+ clc2(e2k5ij  + e2j5ik  " e2i5jk)I2 
+ C2e1ie1jelkdl3  + C2(e2ie1jelk 

+ elie2 jelk  + elielje2k)I4 
+ C2(e2ielje2k  + elie2je2k  + e2ie2jelk)l5 
+ c2e2ie2je2k  I6  > 

where 

= ” ( 6 a)  ell  + (Vb_,?a)  / 

1 2 = r;btan«b  - r,a tan^a  - ln(-b)e1/e  , 

ra 


13  = 2riaCOS29a  - 2r)  bCOS2  6 b 

“ {(^b-e2)  + sin0bsin0a  - cos<?  asin0  a)  eu  , 

14  = 2>?bsin0bcos0b  “ 2»jasin^acos^a 

- (sin20b  - sin20a)eli 

15  = 2r?bsin20b  - 2ry  asin2  9 a - (0b"0a)el  i 

+ (sin0bcos0b  - sinfl acos0 a) eli 

and 

Ig  = 2r?btan0bsin20b  - 2r?atan0  asin2^  a - 2ln(-a)elyg 

, 2 ra 

- ( cos^  9 b - cos^a )elj?; 

O Q 

where  r?a  and  r?b  are  D — evaluated  at  points  a and  b, 

aEe 

respectively . 
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